Prepare data
Merge prevalence data
# create sequence of dates
date_sequence <- seq.Date(min(df_ger_prev$date),
max(df_ger_prev$date), 1)
# create data frame with time sequence
df_dates = tibble(date_sequence, 1:length(date_sequence))
names(df_dates) <- c('date', 'time')
# merge prevalence data
df_ger_prev <- df_ger_prev %>%
inner_join(df_ger_pers, by = 'kreis') %>%
inner_join(df_ger_ctrl, by = 'kreis') %>%
merge(df_dates, by='date') %>%
arrange(kreis, date) %>%
drop_na()
df_ger_prev %>% select(kreis) %>% distinct() %>% nrow()
## [1] 383
# join data frames
df_ger_death <- df_ger_death %>%
plyr::join(df_ger_ctrl, by='kreis') %>%
plyr::join(df_ger_pers, by='kreis') %>%
drop_na()
# save df
df_ger_death %>%
#select(kreis, case_rate, death_rate) %>%
write_csv('ger_death_kreis_controls.csv')
Merge social distancing data
# create sequence of dates
date_sequence <- seq.Date(min(df_ger_socdist$date),
max(df_ger_socdist$date), 1)
# create data frame with time sequence
df_dates = tibble(date_sequence, 1:length(date_sequence))
names(df_dates) <- c('date', 'time')
# merge socdist data
df_ger_socdist <- df_ger_socdist %>%
inner_join(df_ger_pers, by = 'kreis') %>%
inner_join(df_ger_ctrl, by = 'kreis') %>%
merge(df_dates, by='date') %>%
arrange(kreis)
df_ger_socdist %>% select(kreis) %>% distinct() %>% nrow()
## [1] 379
Control for weekend effect
easter <- seq.Date(as.Date('2020-04-10'), as.Date('2020-04-13'), 1)
df_ger_loess <- df_ger_socdist %>%
mutate(weekday = format(date, '%u')) %>%
filter(!(weekday %in% c('6','7') | date %in% easter)) %>%
split(.$kreis) %>%
map(~ loess(socdist_single_tile ~ time, data = .)) %>%
map(predict, 1:max(df_ger_socdist$time)) %>%
bind_rows() %>%
gather(key = 'kreis', value = 'loess') %>%
group_by(kreis) %>%
mutate(time = row_number())
df_ger_loess_2 <- df_ger_socdist %>%
mutate(weekday = format(date, '%u')) %>%
filter(!(weekday %in% c('6','7') | date %in% easter)) %>%
split(.$kreis) %>%
map(~ loess(socdist_tiles ~ time, data = .)) %>%
map(predict, 1:max(df_ger_socdist$time)) %>%
bind_rows() %>%
gather(key = 'kreis', value = 'loess') %>%
rename(loess_2 = loess) %>%
group_by(kreis) %>%
mutate(time = row_number())
df_ger_socdist <- df_ger_socdist %>%
merge(df_ger_loess, by=c('kreis', 'time')) %>%
merge(df_ger_loess_2, by=c('kreis', 'time')) %>%
mutate(weekday = format(date, '%u')) %>%
mutate(socdist_single_tile_clean = ifelse(weekday %in% c('6','7') | date %in% easter,
loess, socdist_single_tile),
socdist_tiles_clean = ifelse(weekday %in% c('6','7') | date %in% easter,
loess_2, socdist_tiles)) %>%
arrange(kreis, time) %>%
select(-weekday)
df_ger_socdist <- df_ger_socdist %>% drop_na() %>% mutate(time = time-1)
Include cleaned social distancing data
df_ger_socdist <- df_ger_socdist %>%
mutate(socdist_single_tile = socdist_single_tile_clean,
socdist_tiles = socdist_tiles_clean) %>%
select(-loess, -loess_2, -socdist_single_tile_clean, -socdist_tiles_clean)
Define outcomes of interest
Social Distancing - Change point analysis
# keep only counties with full data
kreis_complete <- df_ger_socdist %>%
group_by(kreis) %>%
summarise(n = n()) %>%
filter(n==max(.$n)) %>%
.$kreis
# run changepoint analysis
df_ger_socdist_cpt_results <- df_ger_socdist %>%
select(kreis, socdist_single_tile) %>%
filter(kreis %in% kreis_complete) %>%
split(.$kreis) %>%
map(~ cpt.meanvar(as.vector(.$socdist_single_tile),
#penalty = 'Asymptotic',
class=TRUE,
param.estimates=TRUE,
Q=1,
test.stat = 'Normal'))
df_ger_socdist_cpt_results_2 <- df_ger_socdist %>%
select(kreis, socdist_tiles) %>%
filter(kreis %in% kreis_complete) %>%
split(.$kreis) %>%
map(~ cpt.meanvar(as.vector(.$socdist_tiles),
class=TRUE,
param.estimates=TRUE,
Q=1,
test.stat = 'Normal'))
# calculate change point
df_ger_socdist_cpt_day <- df_ger_socdist_cpt_results %>%
map(cpts) %>%
unlist() %>%
as.data.frame() %>%
rename(cpt_day_socdist = '.') %>%
rownames_to_column('kreis')
df_ger_socdist_cpt_day_2 <- df_ger_socdist_cpt_results_2 %>%
map(cpts) %>%
unlist() %>%
as.data.frame() %>%
rename(cpt_day_socdist_2 = '.') %>%
rownames_to_column('kreis')
# calculate mean differences
df_ger_socdist_cpt_mean_diff <- df_ger_socdist_cpt_results %>%
map(param.est) %>%
map(~ .$mean) %>%
map(~ .[2]) %>%
unlist() %>%
as.data.frame() %>%
rename(mean_diff_socdist = '.') %>%
rownames_to_column('kreis')
df_ger_socdist_cpt_mean_diff_2 <- df_ger_socdist_cpt_results_2 %>%
map(param.est) %>%
map(~ .$mean) %>%
map(~ -.[2]) %>%
unlist() %>%
as.data.frame() %>%
rename(mean_diff_socdist_2 = '.') %>%
rownames_to_column('kreis')
# calculate means
df_ger_socdist_mean <- df_ger_socdist %>%
group_by(kreis) %>%
summarise(mean_socdist = mean(socdist_single_tile))
df_ger_socdist_mean_2 <- df_ger_socdist %>%
group_by(kreis) %>%
summarise(mean_socdist_2 = -mean(socdist_tiles))
# merge with county data
df_ger_cpt_socdist <- df_ger_socdist %>%
select(-date, -time, -socdist_single_tile, -socdist_tiles) %>%
distinct() %>%
left_join(df_ger_socdist_cpt_day, by='kreis') %>%
left_join(df_ger_socdist_cpt_day_2, by='kreis') %>%
left_join(df_ger_socdist_cpt_mean_diff, by='kreis') %>%
left_join(df_ger_socdist_cpt_mean_diff_2, by='kreis') %>%
left_join(df_ger_socdist_mean, by='kreis') %>%
left_join(df_ger_socdist_mean_2, by='kreis') %>%
left_join(select(df_ger_onset_prev, kreis, onset_prev), by='kreis') %>%
left_join(select(df_ger_slope_prev, kreis, slope_prev), by='kreis')
# handle censored data
df_ger_cpt_socdist <- df_ger_cpt_socdist %>%
mutate(cpt_day_socdist = ifelse(is.na(cpt_day_socdist),
as.numeric(diff(range(df_ger_socdist$date))),
cpt_day_socdist)) %>%
mutate(event = ifelse(cpt_day_socdist >=
as.numeric(diff(range(df_ger_socdist$date))), 0, 1))
Remove incomplete cases
df_ger_prev_unscaled <- inner_join(df_ger_slope_prev, df_ger_cpt_socdist %>% select(kreis), by='kreis')
df_ger_death_unscaled <- inner_join(df_ger_death, df_ger_cpt_socdist %>% select(kreis), by='kreis')
df_ger_socdist_unscaled <- inner_join(df_ger_cpt_socdist, df_ger_slope_prev %>% select(kreis), by='kreis')
kreis_list <- df_ger_prev_unscaled %>% select(kreis)
kreis_list %>% write_csv('kreis_list.csv')
df_ger_prev_unscaled %>% merge(df_ger_death_unscaled %>% select(kreis:death_rate_0930)) %>%
merge(df_ger_socdist_unscaled %>% select(kreis, cpt_day_socdist_2, mean_diff_socdist_2)) %>%
select(-event) %>% write_csv('ger_all_variables.csv')
Rescale Data
df_ger_prev_scaled <- df_ger_prev_unscaled %>%
mutate_at(vars(-kreis, -event), scale)
df_ger_death_scaled <- df_ger_death_unscaled %>%
mutate_at(vars(-kreis), scale)
df_ger_socdist_scaled <- df_ger_socdist_unscaled %>%
mutate_at(vars(-kreis, -event), scale)
Calculate spatial lags
# create UDF to calculate lagged variables
calc_lags <- function(df, weights, cols){
cols_only <- df %>% select(cols)
cols_lag <- weights %*% as.matrix(cols_only) %>%
as.matrix() %>% as.data.frame()
names(cols_lag) <- paste0(names(cols_lag), '_lag')
return(cols_lag)
}
# read_weight matrix
weight_mat_norm <- read_csv('GER_spatial_weights.csv')
weight_mat_norm <- weight_mat_norm %>% select(-kreis) %>% as.matrix()
dim(weight_mat_norm)
## [1] 379 379
# generate spatially lagged y
y_only <- df_ger_prev_scaled %>% select(onset_prev,slope_prev,slope_prev_var) %>% names()
y_lag <- calc_lags(df_ger_prev_scaled, weight_mat_norm, y_only)
# generate spatially lagged X
X_only <- df_ger_prev_scaled %>% select(academics:pers_n) %>% names()
X_lag <- calc_lags(df_ger_prev_scaled, weight_mat_norm, X_only)
# bind new variables to df
df_ger_prev_scaled <- cbind(df_ger_prev_scaled, y_lag, X_lag)
# generate spatially lagged y
y_only <- df_ger_death_scaled %>% select(case_rate_0331:death_rate_0930) %>% names()
y_lag <- calc_lags(df_ger_death_scaled, weight_mat_norm, y_only)
# generate spatially lagged X
X_only <- df_ger_death_scaled %>% select(academics:pers_n) %>% names()
X_lag <- calc_lags(df_ger_death_scaled, weight_mat_norm, X_only)
# bind new variables to df
df_ger_death_scaled <- cbind(df_ger_death_scaled, y_lag, X_lag)
# generate spatially lagged y
y_only <- df_ger_socdist_scaled %>% select(cpt_day_socdist_2,mean_diff_socdist_2) %>% names()
y_lag <- calc_lags(df_ger_socdist_scaled, weight_mat_norm, y_only)
# generate spatially lagged X
X_only <- df_ger_socdist_scaled %>% select(academics:pers_n) %>% names()
X_lag <- calc_lags(df_ger_socdist_scaled, weight_mat_norm, X_only)
# bind new variables to df
df_ger_socdist_scaled <- cbind(df_ger_socdist_scaled, y_lag, X_lag)
write_csv(df_ger_prev_scaled, 'df_ger_slope_prev.csv')
write_csv(df_ger_death_scaled, 'df_ger_death_scaled.csv')
write_csv(df_ger_socdist_scaled, 'df_ger_cpt_socdist.csv')
Run Specification Curve Analysis
# define function to calculate specification curve analyses for all traits
spec_calculate <- function(df, y, model, controls, all.comb = T){
spec_results_o <- run_specs(df = df,
y = c(y), x = c("pers_o"),
model = c(model), controls = controls,
all.comb = all.comb)
spec_results_c <- run_specs(df = df,
y = c(y), x = c("pers_c"),
model = c(model), controls = controls,
all.comb = all.comb)
spec_results_e <- run_specs(df = df,
y = c(y), x = c("pers_e"),
model = c(model), controls = controls,
all.comb = all.comb)
spec_results_a <- run_specs(df = df,
y = c(y), x = c("pers_a"),
model = c(model), controls = controls,
all.comb = all.comb)
spec_results_n <- run_specs(df = df,
y = c(y), x = c("pers_n"),
model = c(model), controls = controls,
all.comb = all.comb)
spec_results <- list(spec_results_o, spec_results_c, spec_results_e,
spec_results_a, spec_results_n)
names(spec_results) <- list('spec_results_o', 'spec_results_c', 'spec_results_e',
'spec_results_a', 'spec_results_n')
return(spec_results)
}
# adapted from specr package code
format_results <- function(df, var, null = 0, desc = FALSE) {
# rank specs
if (isFALSE(desc)) {
df <- df %>%
dplyr::arrange(!! var)
} else {
df <- df %>%
dplyr::arrange(desc(!! var))
}
# create rank variable and color significance
df <- df %>%
dplyr::mutate(specifications = 1:nrow(df),
color = case_when(conf.low > null ~ "#377eb8",
conf.high < null ~ "#e41a1c",
TRUE ~ "darkgrey"))
return(df)
}
# define function to plot single specification curve
plot_curves <- function(results, filename, hr=F){
file_path_eps <- paste0("../../Plots/GER/", filename,".eps")
file_path_pdf <- paste0("../../Plots/GER/", filename,".pdf")
file_path_png <- paste0("../../Plots/GER/", filename,".png")
if(hr==F){
results %>%
format_results(var = .$estimate, null = 0, desc = F) %>%
ggplot(aes(x = specifications,
y = estimate,
ymin = conf.low,
ymax = conf.high,
color = color)) +
geom_point(aes(color = color),
size = 1) +
theme_minimal() +
scale_color_identity() +
theme(strip.text = element_blank(),
axis.line = element_line("black", size = .5),
legend.position = "none",
panel.spacing = unit(.75, "lines"),
axis.text = element_text(colour = "black")) +
labs(x = "") +
geom_pointrange(alpha = 0.05,
size = .6,
fatten = 1) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
labs(x = "", y = "standarized coefficient") +
coord_fixed(ratio = 2000, ylim = c(-0.4, 0.4)) +
theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black")) +
ggsave(file=file_path_eps, device = 'eps')+
ggsave(file=file_path_pdf, device = 'pdf')+
ggsave(file=file_path_png, device = 'png')
}else{
results %>%
format_results(var = .$estimate, null = 1, desc = F) %>%
ggplot(aes(x = specifications,
y = estimate,
ymin = conf.low,
ymax = conf.high,
color = color)) +
geom_point(aes(color = color),
size = 1) +
theme_minimal() +
scale_color_identity() +
theme(strip.text = element_blank(),
axis.line = element_line("black", size = .5),
legend.position = "none",
panel.spacing = unit(.75, "lines"),
axis.text = element_text(colour = "black")) +
labs(x = "") +
geom_pointrange(alpha = 0.05,
size = .6,
fatten = 1) +
geom_hline(yintercept = 1, linetype = "dashed", color = "black") +
labs(x = "", y = "standarized coefficient") +
coord_fixed(ratio = 2000, ylim = c(0.6, 1.4)) +
theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black")) +
ggsave(file=file_path_eps, device = 'eps')+
ggsave(file=file_path_pdf, device = 'pdf')+
ggsave(file=file_path_png, device = 'png')
}
}
# define function to plot all curves in list
plot_all_curves <- function(ls, analysis, hr=F){
pers <- c('o', 'c', 'e', 'a', 'n')
filenames <- as.list(paste0(analysis, '_', pers))
hr <- as.list(rep(hr, 5))
pmap(list(ls, filenames, hr), plot_curves)
}
# define function to calculate summary stats for single specification curve
calc_summary <- function(df){
dft <- df %>% select(estimate:fit_nobs)
dft <- dft %>% mutate(significant = as.numeric(p.value<0.05),
positive = as.numeric(statistic>0),
negative = as.numeric(statistic<0),
significant_positive = as.numeric(p.value<0.05 & statistic>0),
significant_negative = as.numeric(p.value<0.05 & statistic<0))
mean_temp <- dft %>% map_if(is.numeric, mean, .else=NULL) %>% as.data.frame()
sd_temp <- dft %>% map_if(is.numeric, sd, .else=NULL) %>% as.data.frame()
df_temp <- rbind(mean_temp, sd_temp)
row.names(df_temp) <- c('mean', 'sd')
return(df_temp)
}
# define function to calculate summary stats for all curves in list
calc_all_summaries <- function(ls){
ls <- ls %>% map(calc_summary)
names(ls) <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')
return(ls)
}
coef_to_hr <- function(df){
df %>% mutate(conf.low = exp(estimate-1.96*std.error),
conf.high = exp(estimate+1.96*std.error),
estimate = exp(estimate))
}
filter_socdist <- function(df){
df %>% filter(str_detect(controls, 'onset_prev') &
str_detect(controls, 'slope_prev'))
}
covariates <- c("airport_dist", "conservative", 'male', 'age', 'popdens',
'manufact', 'tourism', 'academics', 'medinc', 'healthcare',
'temp')
Predict Onset
cox_model <- function(formula, data){
formula <- as.formula(formula)
coxph(formula = formula, data = data)}
cox_onset_prev <- spec_calculate(df = df_ger_prev_scaled,
y = "Surv(onset_prev, event)",
model = "cox_model",
controls = covariates %>% append('onset_prev_lag'),
all.comb = T)
cox_onset_prev_hr <- cox_onset_prev %>% map(coef_to_hr)
calc_all_summaries(cox_onset_prev_hr)
## $pers_o
## estimate std.error statistic p.value conf.low conf.high fit_n
## mean 1.02726691 0.055320547 0.4779412 0.5231112 0.92180104 1.14483537 379
## sd 0.05390301 0.002840416 0.9535330 0.2800335 0.05018495 0.05845821 0
## fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean 379 131.39127 3.811112e-06 140.8145
## sd 0 36.26568 5.698238e-05 41.9226
## fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean 4.156490e-06 137.52997 4.234893e-06 NA
## sd 5.770949e-05 40.56007 5.819916e-05 NA
## fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean NA 0.2895930 0.9999496 0.7195157
## sd NA 0.0722395 0.0000000 0.0381976
## fit_std.error.concordance fit_logLik fit_AIC fit_BIC fit_nobs
## mean 0.015635755 -1809.51851 3633.0370 3660.59978 379
## sd 0.001317056 18.13284 34.3554 31.38233 0
## significant positive negative significant_positive significant_negative
## mean 0.08911133 0.6496582 0.3503418 0.08911133 0
## sd 0.28493915 0.4771352 0.4771352 0.28493915 0
##
## $pers_c
## estimate std.error statistic p.value conf.low conf.high fit_n
## mean 0.94922252 0.053595445 -0.9954995 0.3479716 0.85456571 1.05437271 379
## sd 0.04503002 0.001467328 0.8933022 0.2887067 0.04054931 0.05018395 0
## fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean 379 132.01667 0.0001954177 140.8434
## sd 0 37.85635 0.0029957287 43.1195
## fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean 0.0001852582 137.59893 0.0001881611 NA
## sd 0.0028016025 41.78484 0.0028121428 NA
## fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean NA 0.29044877 0.9999496 0.71826746
## sd NA 0.07571748 0.0000000 0.04121107
## fit_std.error.concordance fit_logLik fit_AIC fit_BIC fit_nobs
## mean 0.015744358 -1809.20581 3632.41162 3659.9744 379
## sd 0.001364156 18.92817 35.97484 33.0281 0
## significant positive negative significant_positive significant_negative
## mean 0.1445312 0.1584473 0.8415527 0 0.1445312
## sd 0.3516705 0.3652045 0.3652045 0 0.3516705
##
## $pers_e
## estimate std.error statistic p.value conf.low conf.high fit_n
## mean 1.10484855 0.052014824 1.9131109 0.1418816 0.99784541 1.22333370 379
## sd 0.05117709 0.001260917 0.9230424 0.1475011 0.04805772 0.05449954 0
## fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean 379 134.80506 6.667657e-08 144.44569
## sd 0 33.69771 7.667802e-07 39.58589
## fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean 1.086854e-07 141.00631 1.184551e-07 NA
## sd 1.098564e-06 38.16651 1.180765e-06 NA
## fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean NA 0.29642439 0.9999496 0.72680567
## sd NA 0.06635634 0.0000000 0.03214107
## fit_std.error.concordance fit_logLik fit_AIC fit_BIC fit_nobs
## mean 0.015317233 -1807.81162 3629.62324 3657.18599 379
## sd 0.001142571 16.84885 31.78596 28.87968 0
## significant positive negative significant_positive significant_negative
## mean 0.3879395 1 0 0.3879395 0
## sd 0.4873401 0 0 0.4873401 0
##
## $pers_a
## estimate std.error statistic p.value conf.low conf.high fit_n
## mean 1.006651 0.053837899 0.1292402 0.645567 0.90592390 1.11859211 379
## sd 0.032338 0.001749636 0.6071196 0.238925 0.03165915 0.03294262 0
## fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean 379 130.6175 0.0006586847 140.37322
## sd 0 37.6386 0.0123791555 43.20648
## fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean 0.0006630851 137.15705 0.000668293 NA
## sd 0.0122500086 41.79493 0.012237339 NA
## fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean NA 0.2878630 0.9999496 0.71879483
## sd NA 0.0756924 0.0000000 0.04050933
## fit_std.error.concordance fit_logLik fit_AIC fit_BIC fit_nobs
## mean 0.015616239 -1809.9054 3633.8108 3661.37351 379
## sd 0.001276486 18.8193 35.7401 32.76038 0
## significant positive negative significant_positive significant_negative
## mean 0 0.4870605 0.5129395 0 0
## sd 0 0.4998936 0.4998936 0 0
##
## $pers_n
## estimate std.error statistic p.value conf.low conf.high fit_n
## mean 0.8618104 0.0536108475 -2.7813176 0.01334193 0.77585634 0.95728861 379
## sd 0.0256510 0.0006767025 0.5515346 0.01683804 0.02333785 0.02824371 0
## fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean 379 138.42138 8.297546e-07 148.34185
## sd 0 35.84265 1.236878e-05 41.64845
## fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean 9.623814e-07 144.57863 1.004892e-06 NA
## sd 1.340854e-05 40.36034 1.343165e-05 NA
## fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean NA 0.30271384 0.9999496 0.72815172
## sd NA 0.07048457 0.0000000 0.03346175
## fit_std.error.concordance fit_logLik fit_AIC fit_BIC fit_nobs
## mean 0.015109031 -1806.00346 3626.00692 3653.56967 379
## sd 0.001124059 17.92132 33.95795 31.05319 0
## significant positive negative significant_positive significant_negative
## mean 0.9492188 0 1 0 0.9492188
## sd 0.2195775 0 0 0 0.2195775
plot_all_curves(cox_onset_prev_hr, 'cox_onset_prev_hr', hr = T)
## $spec_results_o
##
## $spec_results_c
##
## $spec_results_e
##
## $spec_results_a
##
## $spec_results_n
cox_onset_prev_ctrl <- coxph(Surv(onset_prev, event) ~ airport_dist +
conservative + male + age + popdens + manufact +
tourism + academics + medinc + healthcare +
temp + onset_prev_lag,
data = df_ger_prev_scaled)
summary(cox_onset_prev_ctrl)
## Call:
## coxph(formula = Surv(onset_prev, event) ~ airport_dist + conservative +
## male + age + popdens + manufact + tourism + academics + medinc +
## healthcare + temp + onset_prev_lag, data = df_ger_prev_scaled)
##
## n= 379, number of events= 379
##
## coef exp(coef) se(coef) z Pr(>|z|)
## airport_dist -0.09147 0.91259 0.07045 -1.298 0.194133
## conservative -0.08051 0.92264 0.08006 -1.006 0.314608
## male -0.13923 0.87003 0.07451 -1.869 0.061664 .
## age -0.25167 0.77750 0.09118 -2.760 0.005779 **
## popdens -0.09156 0.91251 0.07416 -1.235 0.216965
## manufact 0.00861 1.00865 0.05321 0.162 0.871458
## tourism -0.11689 0.88968 0.05988 -1.952 0.050933 .
## academics 0.07279 1.07551 0.08113 0.897 0.369583
## medinc 0.29682 1.34558 0.07898 3.758 0.000171 ***
## healthcare -0.06525 0.93683 0.06938 -0.940 0.347005
## temp -0.07546 0.92731 0.07139 -1.057 0.290488
## onset_prev_lag -0.36945 0.69111 0.08959 -4.124 3.73e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## airport_dist 0.9126 1.0958 0.7949 1.0477
## conservative 0.9226 1.0838 0.7886 1.0794
## male 0.8700 1.1494 0.7518 1.0068
## age 0.7775 1.2862 0.6503 0.9296
## popdens 0.9125 1.0959 0.7891 1.0553
## manufact 1.0086 0.9914 0.9087 1.1195
## tourism 0.8897 1.1240 0.7912 1.0005
## academics 1.0755 0.9298 0.9174 1.2609
## medinc 1.3456 0.7432 1.1526 1.5709
## healthcare 0.9368 1.0674 0.8177 1.0733
## temp 0.9273 1.0784 0.8062 1.0666
## onset_prev_lag 0.6911 1.4469 0.5798 0.8238
##
## Concordance= 0.755 (se = 0.015 )
## Likelihood ratio test= 177.8 on 12 df, p=<2e-16
## Wald test = 183.9 on 12 df, p=<2e-16
## Score (logrank) test = 193.2 on 12 df, p=<2e-16
Predict Slopes
# fixed time windows
lm_slope_prev <- spec_calculate(df = df_ger_prev_scaled,
y = "slope_prev",
model = "lm",
controls = covariates %>% append('slope_prev_lag'),
all.comb = T)
calc_all_summaries(lm_slope_prev)
## $pers_o
## estimate std.error statistic p.value conf.low conf.high
## mean 0.04770701 0.047105679 0.9732636 0.4496175 -0.04492057 0.14033459
## sd 0.04866155 0.003961232 0.9445598 0.3249956 0.04463451 0.05352556
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3762331 0.3647529 0.79419554 34.39651 0.0000103227
## sd 0.1110695 0.1112149 0.06709402 12.97273 0.0001987136
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -445.11260 908.22521 943.66303 235.78389 371.000000
## sd 1.732262 31.40818 60.85522 57.41092 41.98425 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.1411133 0.8693848 0.1306152 0.1411133
## sd 0 0.3481809 0.3370202 0.3370202 0.3481809
## significant_negative
## mean 0
## sd 0
##
## $pers_c
## estimate std.error statistic p.value conf.low conf.high
## mean -0.02472864 0.042382347 -0.5866914 0.4324615 -0.10806835 0.05861107
## sd 0.04143443 0.003375866 0.9654986 0.2623563 0.04195917 0.04196624
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3747063 0.3632077 0.79495961 34.26618 0.001065338
## sd 0.1157000 0.1159113 0.06951753 13.36336 0.018659191
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -445.39270 908.78539 944.22322 236.36101 371.000000
## sd 1.732262 32.35788 62.76145 59.31843 43.73461 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.09521484 0.3103027 0.6896973 0
## sd 0 0.29354730 0.4626740 0.4626740 0
## significant_negative
## mean 0.09521484
## sd 0.29354730
##
## $pers_e
## estimate std.error statistic p.value conf.low conf.high
## mean 0.09475842 0.042916475 2.145207 0.1152644 0.01036841 0.17914843
## sd 0.05415151 0.002851794 1.072831 0.1384433 0.04912212 0.05928462
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3835310 0.3721721 0.78983225 35.33176 6.972021e-08
## sd 0.1038293 0.1038635 0.06319739 12.45035 1.032508e-06
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -443.15352 904.30703 939.74486 233.02529 371.000000
## sd 1.732262 29.85157 57.72422 54.26548 39.24747 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.4655762 1 0 0.4655762
## sd 0 0.4988745 0 0 0.4988745
## significant_negative
## mean 0
## sd 0
##
## $pers_a
## estimate std.error statistic p.value conf.low conf.high
## mean 0.04444294 0.042531589 1.0138143 0.4314829 -0.03919023 0.12807612
## sd 0.03737716 0.003433074 0.8010556 0.3214330 0.03396617 0.04161167
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3757420 0.3642556 0.7944327 34.36378 0.0001083107
## sd 0.1128081 0.1129825 0.0679861 13.13045 0.0019291217
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -445.19538 908.3908 943.82859 235.96953 371.000000
## sd 1.732262 31.75005 61.5464 58.11313 42.64145 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.1437988 0.9528809 0.04711914 0.1437988
## sd 0 0.3509285 0.2119195 0.21191954 0.3509285
## significant_negative
## mean 0
## sd 0
##
## $pers_n
## estimate std.error statistic p.value conf.low conf.high
## mean -0.06378630 0.042452946 -1.458345 0.3027720 -0.14726483 0.01969223
## sd 0.04686406 0.003144902 0.979003 0.2759218 0.05132561 0.04283259
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3785076 0.3670643 0.79281621 34.71477 2.890617e-06
## sd 0.1093210 0.1094544 0.06617581 12.90910 5.052404e-05
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -444.48386 906.96772 942.40555 234.92414 371.000000
## sd 1.732262 31.04716 60.14311 56.72546 41.32336 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.3442383 0 1 0
## sd 0 0.4751772 0 0 0
## significant_negative
## mean 0.3442383
## sd 0.4751772
plot_all_curves(lm_slope_prev, 'lm_slope_prev')
## $spec_results_o
##
## $spec_results_c
##
## $spec_results_e
##
## $spec_results_a
##
## $spec_results_n
lm_slope_prev_ctrl <- lm(slope_prev ~ airport_dist +
conservative + male + age + popdens + manufact +
tourism + academics + medinc + healthcare +
temp + slope_prev_lag,
data = df_ger_prev_scaled)
summary(lm_slope_prev_ctrl)
##
## Call:
## lm(formula = slope_prev ~ airport_dist + conservative + male +
## age + popdens + manufact + tourism + academics + medinc +
## healthcare + temp + slope_prev_lag, data = df_ger_prev_scaled)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.5022 -0.4555 -0.1231 0.3315 3.8222
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.007574 0.036347 -0.208 0.835056
## airport_dist -0.010112 0.046450 -0.218 0.827784
## conservative 0.007307 0.056396 0.130 0.896981
## male -0.085540 0.053273 -1.606 0.109207
## age -0.296574 0.069586 -4.262 2.58e-05 ***
## popdens -0.116212 0.058141 -1.999 0.046368 *
## manufact 0.061335 0.042127 1.456 0.146268
## tourism -0.126881 0.041916 -3.027 0.002645 **
## academics 0.025668 0.056486 0.454 0.649799
## medinc 0.239428 0.052239 4.583 6.29e-06 ***
## healthcare -0.018543 0.047757 -0.388 0.698034
## temp -0.174761 0.049146 -3.556 0.000426 ***
## slope_prev_lag 0.432857 0.060688 7.133 5.31e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7073 on 366 degrees of freedom
## Multiple R-squared: 0.5156, Adjusted R-squared: 0.4997
## F-statistic: 32.47 on 12 and 366 DF, p-value: < 2.2e-16
# variable time windows
lm_slope_prev_var <- spec_calculate(df = df_ger_prev_scaled,
y = "slope_prev_var",
model = "lm",
controls = covariates %>% append('slope_prev_var_lag'),
all.comb = T)
calc_all_summaries(lm_slope_prev_var)
## $pers_o
## estimate std.error statistic p.value conf.low conf.high
## mean 0.03138877 0.048027782 0.6168286 0.5240564 -0.06305202 0.1258296
## sd 0.04551793 0.003904359 0.8881225 0.3037374 0.04165726 0.0502624
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3513107 0.3393623 0.81020283 30.71510 0.0001379822
## sd 0.1089591 0.1089781 0.06488532 11.75029 0.0027232807
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -452.79671 923.59341 959.03124 245.20455 371.000000
## sd 1.732262 29.99134 57.97282 54.44501 41.18655 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.07739258 0.7424316 0.2575684 0.07739258
## sd 0 0.26724596 0.4373484 0.4373484 0.26724596
## significant_negative
## mean 0
## sd 0
##
## $pers_c
## estimate std.error statistic p.value conf.low conf.high
## mean -0.005412083 0.043251658 -0.1122928 0.4938687 -0.09046118 0.07963702
## sd 0.041568173 0.003203958 0.9510700 0.2765537 0.04319003 0.04086338
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3505259 0.3385690 0.81057266 30.65762 0.001475911
## sd 0.1118288 0.1118893 0.06636315 11.99597 0.025347073
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -452.92184 923.8437 959.28150 245.50122 371.000000
## sd 1.732262 30.55817 59.1121 55.58766 42.27129 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.04272461 0.4780273 0.5219727 0
## sd 0 0.20226024 0.4995780 0.4995780 0
## significant_negative
## mean 0.04272461
## sd 0.20226024
##
## $pers_e
## estimate std.error statistic p.value conf.low conf.high
## mean 0.06672988 0.043867749 1.464399 0.3014781 -0.01953070 0.1529904
## sd 0.05252357 0.002821869 1.068661 0.2892869 0.04787854 0.0573294
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3555855 0.3437053 0.80774108 31.20591 2.225108e-06
## sd 0.1036970 0.1036339 0.06204795 11.33341 3.180183e-05
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -451.73375 921.46749 956.90532 243.58870 371.000000
## sd 1.732262 28.85099 55.67972 52.14376 39.19745 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.2744141 0.9606934 0.03930664 0.2744141
## sd 0 0.4462730 0.1943472 0.19434724 0.4462730
## significant_negative
## mean 0
## sd 0
##
## $pers_a
## estimate std.error statistic p.value conf.low conf.high
## mean 0.05484896 0.043310646 1.2488281 0.3161042 -0.03031614 0.14001406
## sd 0.03452171 0.003320092 0.7388706 0.2662508 0.03219300 0.03784622
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3528556 0.3409354 0.8092101 30.95031 0.0001543231
## sd 0.1093621 0.1094004 0.0651510 11.88643 0.0025729140
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -452.3210 922.64193 958.07975 244.62057 371.000000
## sd 1.732262 30.1185 58.24095 54.74094 41.33889 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.1984863 0.9831543 0.0168457 0.1984863
## sd 0 0.3989090 0.1287089 0.1287089 0.3989090
## significant_negative
## mean 0
## sd 0
##
## $pers_n
## estimate std.error statistic p.value conf.low conf.high
## mean -0.01945947 0.043371756 -0.4091114 0.4680019 -0.10474473 0.06582579
## sd 0.04490624 0.003125042 0.9702821 0.2636699 0.04887248 0.04147441
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3511661 0.3392152 0.81027906 30.70800 7.639954e-05
## sd 0.1092723 0.1093028 0.06506678 11.78065 1.248112e-03
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -452.8261 923.65224 959.09007 245.25920 371.000000
## sd 1.732262 30.0687 58.13108 54.60971 41.30493 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.08203125 0.4973145 0.5026855 0
## sd 0 0.27444583 0.5000538 0.5000538 0
## significant_negative
## mean 0.08203125
## sd 0.27444583
plot_all_curves(lm_slope_prev_var, 'lm_slope_prev_var')
## $spec_results_o
##
## $spec_results_c
##
## $spec_results_e
##
## $spec_results_a
##
## $spec_results_n
lm_slope_prev_var_ctrl <- lm(slope_prev_var ~ airport_dist +
conservative + male + age + popdens + manufact +
tourism + academics + medinc + healthcare +
temp + slope_prev_var_lag,
data = df_ger_prev_scaled)
summary(lm_slope_prev_var_ctrl)
##
## Call:
## lm(formula = slope_prev_var ~ airport_dist + conservative + male +
## age + popdens + manufact + tourism + academics + medinc +
## healthcare + temp + slope_prev_var_lag, data = df_ger_prev_scaled)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.4605 -0.4684 -0.0894 0.3536 3.8629
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0107049 0.0370087 -0.289 0.772552
## airport_dist -0.0006466 0.0472676 -0.014 0.989094
## conservative 0.0578272 0.0574076 1.007 0.314452
## male -0.1063530 0.0542440 -1.961 0.050679 .
## age -0.3348173 0.0704594 -4.752 2.90e-06 ***
## popdens -0.1144899 0.0591642 -1.935 0.053745 .
## manufact 0.0844884 0.0428457 1.972 0.049371 *
## tourism -0.1367335 0.0426814 -3.204 0.001476 **
## academics -0.0251356 0.0576842 -0.436 0.663279
## medinc 0.2276969 0.0521603 4.365 1.65e-05 ***
## healthcare -0.0223008 0.0487212 -0.458 0.647423
## temp -0.1940683 0.0505286 -3.841 0.000145 ***
## slope_prev_var_lag 0.4329906 0.0602916 7.182 3.88e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7199 on 366 degrees of freedom
## Multiple R-squared: 0.4982, Adjusted R-squared: 0.4817
## F-statistic: 30.28 on 12 and 366 DF, p-value: < 2.2e-16
Predict Cases
lm_cases_0331 <- spec_calculate(df = df_ger_death_scaled,
y = "case_rate_0331",
model = "lm",
controls = covariates %>% append('case_rate_0331_lag'),
all.comb = T)
calc_all_summaries(lm_cases_0331)
## $pers_o
## estimate std.error statistic p.value conf.low conf.high
## mean 0.03158699 0.051022780 0.5994671 0.5627758 -0.06874312 0.13191710
## sd 0.04075752 0.003039436 0.7655479 0.2958935 0.03824516 0.04394422
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.27125522 0.2577158 0.86029794 20.714339 0.0002426819
## sd 0.08294606 0.0823706 0.04660655 7.515557 0.0038974652
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -476.16370 970.32739 1005.76522 275.46553 371.000000
## sd 1.732262 20.54957 39.18622 36.10197 31.35361 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.06103516 0.7888184 0.2111816 0.06103516
## sd 0 0.23942402 0.4081968 0.4081968 0.23942402
## significant_negative
## mean 0
## sd 0
##
## $pers_c
## estimate std.error statistic p.value conf.low conf.high
## mean -0.004091242 0.045901847 -0.08596589 0.5796038 -0.09435164 0.08616916
## sd 0.034700872 0.002269187 0.74910889 0.2613804 0.03536478 0.03460409
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.27028449 0.25673269 0.8607943 20.626725 0.002355802
## sd 0.08560184 0.08506666 0.0479699 7.724115 0.031524022
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -476.35290 970.70581 1006.14364 275.8325 371.000000
## sd 1.732262 21.07044 40.23365 37.14286 32.3575 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0 0.4645996 0.5354004 0
## sd 0 0 0.4988061 0.4988061 0
## significant_negative
## mean 0
## sd 0
##
## $pers_e
## estimate std.error statistic p.value conf.low conf.high
## mean 0.08714459 0.046581747 1.8422270 0.1520497 -0.004452756 0.17874193
## sd 0.04647008 0.001763032 0.9112771 0.1493776 0.043558178 0.04945351
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.27812366 0.26469914 0.85640497 21.418135 2.574047e-06
## sd 0.07660416 0.07592723 0.04326494 7.115188 3.140562e-05
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -474.51363 967.0273 1002.46509 272.86926 371.000000
## sd 1.732262 19.24243 36.5560 33.48926 28.95637 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.3618164 1 0 0.3618164
## sd 0 0.4805847 0 0 0.4805847
## significant_negative
## mean 0
## sd 0
##
## $pers_a
## estimate std.error statistic p.value conf.low conf.high
## mean 0.07275316 0.045946320 1.5732114 0.1733925 -0.01759469 0.16310101
## sd 0.02970329 0.002372876 0.5983738 0.1394271 0.02789450 0.03209363
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.27494180 0.26146891 0.85812315 21.120691 0.0001372159
## sd 0.08246331 0.08190365 0.04643572 7.589597 0.0017659528
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -475.20596 968.41191 1003.8497 274.07200 371.000000
## sd 1.732262 20.51326 39.13057 36.0854 31.17113 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.2634277 1 0 0.2634277
## sd 0 0.4405462 0 0 0.4405462
## significant_negative
## mean 0
## sd 0
##
## $pers_n
## estimate std.error statistic p.value conf.low conf.high
## mean -0.06299523 0.045987096 -1.352350 0.2837433 -0.15342327 0.02743280
## sd 0.03825584 0.002107881 0.770264 0.2261862 0.04073718 0.03608117
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.27430368 0.26081464 0.85855349 21.04466 2.747966e-05
## sd 0.08078661 0.08020112 0.04551672 7.43662 3.526301e-04
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -475.41660 968.83321 1004.27103 274.31321 371.000000
## sd 1.732262 20.13405 38.37129 35.33896 30.53734 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.2299805 0 1 0
## sd 0 0.4208714 0 0 0
## significant_negative
## mean 0.2299805
## sd 0.4208714
plot_all_curves(lm_cases_0331, 'lm_cases_0331')
## $spec_results_o
##
## $spec_results_c
##
## $spec_results_e
##
## $spec_results_a
##
## $spec_results_n
lm_cases_0331_ctrl <- lm(case_rate_0331 ~ airport_dist +
conservative + male + age + popdens + manufact +
tourism + academics + medinc + healthcare +
temp + case_rate_0331_lag,
data = df_ger_death_scaled)
summary(lm_cases_0331_ctrl)
##
## Call:
## lm(formula = case_rate_0331 ~ airport_dist + conservative + male +
## age + popdens + manufact + tourism + academics + medinc +
## healthcare + temp + case_rate_0331_lag, data = df_ger_death_scaled)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.4908 -0.4246 -0.1322 0.1935 5.7140
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.005574 0.041192 -0.135 0.892434
## airport_dist 0.007239 0.052639 0.138 0.890700
## conservative 0.021561 0.064001 0.337 0.736398
## male -0.079280 0.060385 -1.313 0.190032
## age -0.287805 0.078147 -3.683 0.000265 ***
## popdens -0.094136 0.065815 -1.430 0.153484
## manufact 0.045203 0.047719 0.947 0.344117
## tourism -0.094058 0.047518 -1.979 0.048517 *
## academics -0.011510 0.063625 -0.181 0.856548
## medinc 0.184173 0.059260 3.108 0.002032 **
## healthcare -0.036733 0.054269 -0.677 0.498913
## temp -0.180424 0.055771 -3.235 0.001327 **
## case_rate_0331_lag 0.404243 0.071769 5.633 3.54e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8017 on 366 degrees of freedom
## Multiple R-squared: 0.3777, Adjusted R-squared: 0.3573
## F-statistic: 18.51 on 12 and 366 DF, p-value: < 2.2e-16
lm_cases_0930 <- spec_calculate(df = df_ger_death_scaled,
y = "case_rate_0930",
model = "lm",
controls = covariates %>% append('case_rate_0930_lag'),
all.comb = T)
calc_all_summaries(lm_cases_0930)
## $pers_o
## estimate std.error statistic p.value conf.low conf.high
## mean 0.04703511 0.048759556 0.9328083 0.4213025 -0.04884462 0.14291485
## sd 0.04623303 0.003708065 0.8870397 0.3049813 0.04263197 0.05063419
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3311692 0.3188131 0.82304322 28.08371 0.0000322385
## sd 0.1042954 0.1042904 0.06154394 11.21308 0.0004799432
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -458.89825 935.79650 971.23433 252.81805 371.000000
## sd 1.732262 28.15626 54.43081 51.21137 39.42368 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.1430664 0.8413086 0.1586914 0.1430664
## sd 0 0.3501833 0.3654327 0.3654327 0.3501833
## significant_negative
## mean 0
## sd 0
##
## $pers_c
## estimate std.error statistic p.value conf.low conf.high
## mean -0.002610167 0.044018978 -0.02673773 0.5244260 -0.08916812 0.08394778
## sd 0.040295483 0.003040096 0.89565674 0.2664612 0.04373607 0.03749754
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3292839 0.3169004 0.82408293 27.87911 0.0006614882
## sd 0.1074214 0.1074516 0.06314984 11.42535 0.0139028188
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -459.32868 936.65736 972.09518 253.53069 371.000000
## sd 1.732262 28.77586 55.66624 52.42724 40.60527 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.0324707 0.5195312 0.4804688 0
## sd 0 0.1772682 0.4996794 0.4996794 0
## significant_negative
## mean 0.0324707
## sd 0.1772682
##
## $pers_e
## estimate std.error statistic p.value conf.low conf.high
## mean 0.11873043 0.044310503 2.6373935 0.03352240 0.03159922 0.20586163
## sd 0.04761501 0.002634405 0.9060562 0.04271978 0.04339097 0.05201209
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.34275305 0.33059563 0.81616073 29.48154 1.913204e-08
## sd 0.09612161 0.09600844 0.05733096 10.80161 3.024722e-07
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -455.83276 929.6655 965.10334 248.43935 371.000000
## sd 1.732262 26.55416 51.2221 48.02528 36.33397 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.7509766 1 0 0.7509766
## sd 0 0.4325002 0 0 0.4325002
## significant_negative
## mean 0
## sd 0
##
## $pers_a
## estimate std.error statistic p.value conf.low conf.high
## mean 0.001807263 0.044142706 0.02930601 0.6298433 -0.08499398 0.08860851
## sd 0.029309078 0.003273494 0.63816227 0.2379830 0.02826462 0.03165452
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.328552 0.3161559 0.8245256 27.78624 0.0008432176
## sd 0.107648 0.1076773 0.0632661 11.42428 0.0168326235
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -459.52910 937.0582 972.49602 253.80733 371.000000
## sd 1.732262 28.81963 55.7506 52.50393 40.69093 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0 0.5290527 0.4709473 0
## sd 0 0 0.4992162 0.4992162 0
## significant_negative
## mean 0
## sd 0
##
## $pers_n
## estimate std.error statistic p.value conf.low conf.high
## mean -0.001840958 0.044147329 -0.01261312 0.5589962 -0.08865129 0.08496938
## sd 0.037168325 0.003113077 0.79081236 0.2465463 0.04075222 0.03430969
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3290563 0.3166673 0.82425385 27.83814 0.0001595196
## sd 0.1066349 0.1066540 0.06276337 11.33910 0.0028612546
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -459.41897 936.83794 972.27576 253.61671 371.000000
## sd 1.732262 28.63422 55.37859 52.13274 40.30798 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.01220703 0.5891113 0.4108887 0
## sd 0 0.10982242 0.4920552 0.4920552 0
## significant_negative
## mean 0.01220703
## sd 0.10982242
plot_all_curves(lm_cases_0930, 'lm_cases_0930')
## $spec_results_o
##
## $spec_results_c
##
## $spec_results_e
##
## $spec_results_a
##
## $spec_results_n
lm_cases_0930_ctrl <- lm(case_rate_0930 ~ airport_dist +
conservative + male + age + popdens + manufact +
tourism + academics + medinc + healthcare +
temp + case_rate_0930_lag,
data = df_ger_death_scaled)
summary(lm_cases_0930_ctrl)
##
## Call:
## lm(formula = case_rate_0930 ~ airport_dist + conservative + male +
## age + popdens + manufact + tourism + academics + medinc +
## healthcare + temp + case_rate_0930_lag, data = df_ger_death_scaled)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6305 -0.4216 -0.1215 0.3040 3.5774
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.003935 0.038008 -0.104 0.917604
## airport_dist -0.045205 0.048579 -0.931 0.352702
## conservative 0.102075 0.058827 1.735 0.083551 .
## male -0.034565 0.055737 -0.620 0.535548
## age -0.350952 0.071996 -4.875 1.63e-06 ***
## popdens 0.152294 0.061336 2.483 0.013478 *
## manufact 0.062070 0.044049 1.409 0.159651
## tourism -0.109885 0.043906 -2.503 0.012760 *
## academics -0.217018 0.059681 -3.636 0.000316 ***
## medinc 0.190831 0.056091 3.402 0.000742 ***
## healthcare 0.013822 0.049981 0.277 0.782285
## temp -0.156873 0.051285 -3.059 0.002385 **
## case_rate_0930_lag 0.479824 0.068459 7.009 1.16e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7399 on 366 degrees of freedom
## Multiple R-squared: 0.47, Adjusted R-squared: 0.4526
## F-statistic: 27.05 on 12 and 366 DF, p-value: < 2.2e-16
Predict Deaths
lm_deaths_0331 <- spec_calculate(df = df_ger_death_scaled,
y = "death_rate_0331",
model = "lm",
controls = covariates %>% append('death_rate_0331_lag'),
all.comb = T)
calc_all_summaries(lm_deaths_0331)
## $pers_o
## estimate std.error statistic p.value conf.low conf.high
## mean 0.05926579 0.056344799 1.048407 0.3391156 -0.05152945 0.17006104
## sd 0.02649895 0.002223663 0.465681 0.2207861 0.02562865 0.02803226
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.11230093 0.09564048 0.95084691 6.817236 0.003741489
## sd 0.03249807 0.03024746 0.01580279 1.831876 0.022251556
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -514.578479 1047.15696 1082.59478 335.55025 371.000000
## sd 1.732262 6.852099 11.52894 9.74724 12.28427 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.02636719 0.9821777 0.01782227 0.02636719
## sd 0 0.16024428 0.1323212 0.13232123 0.16024428
## significant_negative
## mean 0
## sd 0
##
## $pers_c
## estimate std.error statistic p.value conf.low conf.high
## mean -0.02176714 0.0507656910 -0.4272184 0.6489583 -0.12159172 0.07805745
## sd 0.02201730 0.0008507798 0.4327864 0.2340294 0.02246288 0.02169193
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.11002111 0.09331891 0.95206101 6.651521 0.00963241
## sd 0.03320389 0.03098238 0.01615476 1.883209 0.05665206
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -515.059950 1048.11990 1083.557726 336.41202 371.000000
## sd 1.732262 6.973763 11.79201 9.998775 12.55107 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0 0.1816406 0.8183594 0
## sd 0 0 0.3855951 0.3855951 0
## significant_negative
## mean 0
## sd 0
##
## $pers_e
## estimate std.error statistic p.value conf.low conf.high
## mean 0.06219816 0.0515241068 1.2102704 0.2768198 -0.03911776 0.16351408
## sd 0.02481655 0.0007406798 0.4899089 0.1810235 0.02544552 0.02425882
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.11323239 0.09658630 0.95035178 6.898352 0.002436617
## sd 0.03213563 0.02994777 0.01566049 1.848994 0.014602120
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -514.381924 1046.76385 1082.201674 335.19815 371.000000
## sd 1.732262 6.787753 11.45944 9.857845 12.14727 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.09594727 1 0 0.09594727
## sd 0 0.29455487 0 0 0.29455487
## significant_negative
## mean 0
## sd 0
##
## $pers_a
## estimate std.error statistic p.value conf.low conf.high
## mean 0.07596961 0.050837168 1.4956105 0.1615796 -0.02399553 0.17593475
## sd 0.01854741 0.001137139 0.3694743 0.1092908 0.01880931 0.01855325
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.11477591 0.09816324 0.9495095 6.996946 0.003652699
## sd 0.03348276 0.03131897 0.0163866 1.927123 0.022385471
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -514.040955 1046.08191 1081.51974 334.61471 371.000000
## sd 1.732262 7.080364 12.03403 10.26139 12.65648 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.1005859 1 0 0.1005859
## sd 0 0.3008164 0 0 0.3008164
## significant_negative
## mean 0
## sd 0
##
## $pers_n
## estimate std.error statistic p.value conf.low conf.high
## mean -0.05389933 0.050903492 -1.0622345 0.3222672 -0.15399488 0.04619623
## sd 0.01914286 0.001035423 0.3829244 0.1659477 0.01845582 0.02001427
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.11220076 0.09553655 0.95090074 6.820286 0.003463388
## sd 0.03252479 0.03033547 0.01585213 1.860572 0.020486709
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -514.599597 1047.19919 1082.637020 335.58811 371.000000
## sd 1.732262 6.860555 11.59847 9.943315 12.29437 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.005859375 0 1 0
## sd 0 0.076331286 0 0 0
## significant_negative
## mean 0.005859375
## sd 0.076331286
plot_all_curves(lm_deaths_0331, 'lm_deaths_0331')
## $spec_results_o
##
## $spec_results_c
##
## $spec_results_e
##
## $spec_results_a
##
## $spec_results_n
lm_deaths_0331_ctrl <- lm(death_rate_0331 ~ airport_dist +
conservative + male + age + popdens + manufact +
tourism + academics + medinc + healthcare +
temp + death_rate_0331_lag,
data = df_ger_death_scaled)
summary(lm_deaths_0331_ctrl)
##
## Call:
## lm(formula = death_rate_0331 ~ airport_dist + conservative +
## male + age + popdens + manufact + tourism + academics + medinc +
## healthcare + temp + death_rate_0331_lag, data = df_ger_death_scaled)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.8867 -0.4742 -0.2131 0.1644 5.1407
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.003276 0.047481 -0.069 0.945030
## airport_dist 0.063879 0.060650 1.053 0.292924
## conservative 0.033146 0.073530 0.451 0.652415
## male -0.170458 0.069611 -2.449 0.014805 *
## age -0.282595 0.087552 -3.228 0.001360 **
## popdens -0.056660 0.075956 -0.746 0.456170
## manufact 0.190338 0.055019 3.460 0.000605 ***
## tourism -0.071972 0.054777 -1.314 0.189699
## academics -0.081969 0.072606 -1.129 0.259654
## medinc 0.061189 0.065733 0.931 0.352534
## healthcare -0.095838 0.062240 -1.540 0.124469
## temp -0.120116 0.064268 -1.869 0.062423 .
## death_rate_0331_lag 0.231548 0.087866 2.635 0.008766 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.924 on 366 degrees of freedom
## Multiple R-squared: 0.1733, Adjusted R-squared: 0.1462
## F-statistic: 6.392 on 12 and 366 DF, p-value: 2.503e-10
lm_deaths_0930 <- spec_calculate(df = df_ger_death_scaled,
y = "death_rate_0930",
model = "lm",
controls = covariates %>% append('death_rate_0930_lag'),
all.comb = T)
calc_all_summaries(lm_deaths_0930)
## $pers_o
## estimate std.error statistic p.value conf.low conf.high
## mean 0.09176990 0.055963743 1.6313826 0.145225 -0.01827605 0.20181584
## sd 0.02910357 0.002231346 0.4896016 0.137522 0.02703758 0.03164673
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.12375372 0.1073184 0.94465296 7.626281 0.002057778
## sd 0.03576864 0.0336473 0.01767784 2.103161 0.013217159
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -512.087716 1042.17543 1077.6133 331.22109 371.000000
## sd 1.732262 7.626127 13.09662 11.0166 13.52055 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.2592773 1 0 0.2592773
## sd 0 0.4382916 0 0 0.4382916
## significant_negative
## mean 0
## sd 0
##
## $pers_c
## estimate std.error statistic p.value conf.low conf.high
## mean -0.02373502 0.0505987501 -0.4672378 0.6381103 -0.12323134 0.07576129
## sd 0.02148019 0.0009822833 0.4226164 0.2372276 0.02203167 0.02109180
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.11782829 0.10128121 0.94783532 7.197791 0.009895088
## sd 0.03662981 0.03456711 0.01808525 2.168291 0.058343172
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -513.359819 1044.71964 1080.15746 333.46091 371.000000
## sd 1.732262 7.745196 13.38859 11.39777 13.84607 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0 0.1352539 0.8647461 0
## sd 0 0 0.3420363 0.3420363 0
## significant_negative
## mean 0
## sd 0
##
## $pers_e
## estimate std.error statistic p.value conf.low conf.high
## mean 0.05256399 0.051225663 1.0284875 0.3487208 -0.04816508 0.15329305
## sd 0.02270246 0.000942625 0.4474427 0.1962162 0.02322695 0.02232001
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.11989229 0.10338102 0.94673451 7.358761 0.003794686
## sd 0.03579107 0.03374005 0.01768661 2.135709 0.021709173
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -512.922157 1043.8443 1079.28214 332.68071 371.000000
## sd 1.732262 7.596441 13.1105 11.22513 13.52902 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.04541016 1 0 0.04541016
## sd 0 0.20822742 0 0 0.20822742
## significant_negative
## mean 0
## sd 0
##
## $pers_a
## estimate std.error statistic p.value conf.low conf.high
## mean 0.01004584 0.05075540 0.1959914 0.7221615 -0.08975852 0.10985019
## sd 0.02082310 0.00116779 0.4083976 0.1840059 0.02043765 0.02144881
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.11736593 0.1008111 0.94808125 7.161593 0.01013347
## sd 0.03683933 0.0347719 0.01819243 2.173253 0.06061994
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -513.457345 1044.91469 1080.3525 333.63568 371.000000
## sd 1.732262 7.789187 13.46958 11.4449 13.92527 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0 0.6926270 0.3073730 0
## sd 0 0 0.4614616 0.4614616 0
## significant_negative
## mean 0
## sd 0
##
## $pers_n
## estimate std.error statistic p.value conf.low conf.high
## mean 0.03733565 0.050718234 0.7335759 0.4914366 -0.06239562 0.13706691
## sd 0.02251156 0.001073727 0.4357044 0.2497344 0.02195929 0.02324317
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.11858409 0.10205446 0.94741952 7.243279 0.01153319
## sd 0.03739883 0.03533075 0.01848994 2.202537 0.06886308
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -513.190383 1044.38077 1079.81859 333.17521 371.000000
## sd 1.732262 7.912561 13.70082 11.59235 14.13676 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0 0.9448242 0.05517578 0
## sd 0 0 0.2283509 0.22835092 0
## significant_negative
## mean 0
## sd 0
plot_all_curves(lm_deaths_0930, 'lm_deaths_0930')
## $spec_results_o
##
## $spec_results_c
##
## $spec_results_e
##
## $spec_results_a
##
## $spec_results_n
lm_deaths_0930_ctrl <- lm(death_rate_0930 ~ airport_dist +
conservative + male + age + popdens + manufact +
tourism + academics + medinc + healthcare +
temp + death_rate_0930_lag,
data = df_ger_death_scaled)
summary(lm_deaths_0930_ctrl)
##
## Call:
## lm(formula = death_rate_0930 ~ airport_dist + conservative +
## male + age + popdens + manufact + tourism + academics + medinc +
## healthcare + temp + death_rate_0930_lag, data = df_ger_death_scaled)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.9553 -0.5092 -0.2248 0.2556 6.5021
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.004634 0.047233 -0.098 0.92190
## airport_dist 0.026542 0.060328 0.440 0.66023
## conservative 0.026796 0.073079 0.367 0.71407
## male -0.087034 0.069386 -1.254 0.21052
## age -0.211738 0.086909 -2.436 0.01531 *
## popdens 0.029944 0.075712 0.396 0.69270
## manufact 0.135993 0.054799 2.482 0.01353 *
## tourism -0.093115 0.054578 -1.706 0.08884 .
## academics -0.162248 0.072481 -2.238 0.02579 *
## medinc 0.093815 0.065408 1.434 0.15234
## healthcare -0.026707 0.061800 -0.432 0.66588
## temp -0.184670 0.064375 -2.869 0.00436 **
## death_rate_0930_lag 0.284429 0.087564 3.248 0.00127 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9191 on 366 degrees of freedom
## Multiple R-squared: 0.182, Adjusted R-squared: 0.1552
## F-statistic: 6.788 on 12 and 366 DF, p-value: 4.504e-11
Predict Socdist Onset
cox_onset_socdist <- spec_calculate(df = df_ger_socdist_scaled,
y = "Surv(cpt_day_socdist_2, event)",
model = "cox_model",
controls = covariates %>%
append(c('cpt_day_socdist_2_lag', 'onset_prev', 'slope_prev')),
all.comb = T)
cox_onset_socdist <- cox_onset_socdist %>% map(filter_socdist)
cox_onset_socdist_hr <- cox_onset_socdist %>% map(coef_to_hr)
calc_all_summaries(cox_onset_socdist_hr)
## $pers_o
## estimate std.error statistic p.value conf.low conf.high fit_n
## mean 0.92412115 0.060146732 -1.3205350 0.2201715 0.82135075 1.03977081 379
## sd 0.02156262 0.002248289 0.4019257 0.1387493 0.01899451 0.02532624 0
## fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean 379 38.32341 0.005706523 36.55378
## sd 0 13.75696 0.024395795 12.35949
## fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean 0.005739844 36.14392 0.005735019 NA
## sd 0.023985584 12.05654 0.023777053 NA
## fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean NA 0.09557573 0.9999496 0.61749673
## sd NA 0.03296967 0.0000000 0.03336312
## fit_std.error.concordance fit_logLik fit_AIC fit_BIC fit_nobs
## mean 0.028868295 -1856.052442 3730.10488 3765.543 379
## sd 0.001210389 6.878482 12.02107 11.250 0
## significant positive negative significant_positive significant_negative
## mean 0.06811523 0 1 0 0.06811523
## sd 0.25197430 0 0 0 0.25197430
##
## $pers_c
## estimate std.error statistic p.value conf.low conf.high fit_n
## mean 0.92379825 0.053425612 -1.4965299 0.1765796 0.8319259 1.02582137 379
## sd 0.02216087 0.001129165 0.4722299 0.1393372 0.0186606 0.02630592 0
## fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean 379 38.89125 0.002080154 37.09716
## sd 0 12.54155 0.007141225 11.33829
## fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean 0.002219141 36.66871 0.002244855 NA
## sd 0.007425369 11.03882 0.007438921 NA
## fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean NA 0.09703073 0.9999496 0.62825725
## sd NA 0.02998817 0.0000000 0.02806221
## fit_std.error.concordance fit_logLik fit_AIC fit_BIC fit_nobs
## mean 0.02888723 -1855.768524 3729.53705 3764.97487 379
## sd 0.00105630 6.270773 10.78844 10.29674 0
## significant positive negative significant_positive significant_negative
## mean 0.1914062 0 1 0 0.1914062
## sd 0.3934561 0 0 0 0.3934561
##
## $pers_e
## estimate std.error statistic p.value conf.low conf.high fit_n
## mean 0.88750828 0.0562011720 -2.1242714 0.04010842 0.79494811 0.99084875 379
## sd 0.01407696 0.0008904741 0.2702866 0.02493297 0.01325974 0.01506055 0
## fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean 379 41.01862 0.002206119 39.25972
## sd 0 13.53017 0.009699107 12.25642
## fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean 0.00234413 38.82729 0.002322961 NA
## sd 0.01021985 11.92044 0.010047201 NA
## fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean NA 0.1020040 0.9999496 0.62851368
## sd NA 0.0321979 0.0000000 0.02755481
## fit_std.error.concordance fit_logLik fit_AIC fit_BIC fit_nobs
## mean 0.028252338 -1854.704839 3727.40968 3762.84750 379
## sd 0.001407939 6.765084 11.72497 10.85692 0
## significant positive negative significant_positive significant_negative
## mean 0.6994629 0 1 0 0.6994629
## sd 0.4585476 0 0 0 0.4585476
##
## $pers_a
## estimate std.error statistic p.value conf.low conf.high fit_n
## mean 0.94998162 0.055857806 -0.9252410 0.3798245 0.85145243 1.0599184 379
## sd 0.01686576 0.001206494 0.3290921 0.1672532 0.01424408 0.0201586 0
## fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean 379 37.40628 0.00487520 35.73613
## sd 0 13.19216 0.01728311 11.95547
## fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean 0.005074892 35.31908 0.005112128 NA
## sd 0.017780079 11.64557 0.017763039 NA
## fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean NA 0.09343289 0.9999496 0.62513121
## sd NA 0.03167007 0.0000000 0.03199302
## fit_std.error.concordance fit_logLik fit_AIC fit_BIC fit_nobs
## mean 0.028872705 -1856.51101 3731.02201 3766.45984 379
## sd 0.001126387 6.59608 11.46145 10.84048 0
## significant positive negative significant_positive significant_negative
## mean 0 0 1 0 0
## sd 0 0 0 0 0
##
## $pers_n
## estimate std.error statistic p.value conf.low conf.high fit_n
## mean 1.08916744 0.057407948 1.4897329 0.14528486 0.97327553 1.21886580 379
## sd 0.01220407 0.001200459 0.2138959 0.06005665 0.01235991 0.01223976 0
## fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean 379 38.68176 0.00353503 37.02830
## sd 0 13.43207 0.01343931 12.20186
## fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean 0.003639324 36.63856 0.003605216 NA
## sd 0.013488312 11.86728 0.013282388 NA
## fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean NA 0.09645849 0.9999496 0.62898123
## sd NA 0.03214506 0.0000000 0.03045481
## fit_std.error.concordance fit_logLik fit_AIC fit_BIC fit_nobs
## mean 0.028670041 -1855.873268 3729.74654 3765.18436 379
## sd 0.001247607 6.716033 11.72373 11.09026 0
## significant positive negative significant_positive significant_negative
## mean 0 1 0 0 0
## sd 0 0 0 0 0
plot_all_curves(cox_onset_socdist_hr, 'cox_onset_socdist_hr', hr = T)
## $spec_results_o
##
## $spec_results_c
##
## $spec_results_e
##
## $spec_results_a
##
## $spec_results_n
cox_onset_socdist_ctrl <- coxph(Surv(cpt_day_socdist_2, event) ~ airport_dist +
conservative + male + age + popdens + manufact +
tourism + academics + medinc + healthcare +
temp + onset_prev + slope_prev +
cpt_day_socdist_2_lag,
data = df_ger_socdist_scaled)
summary(cox_onset_socdist_ctrl)
## Call:
## coxph(formula = Surv(cpt_day_socdist_2, event) ~ airport_dist +
## conservative + male + age + popdens + manufact + tourism +
## academics + medinc + healthcare + temp + onset_prev + slope_prev +
## cpt_day_socdist_2_lag, data = df_ger_socdist_scaled)
##
## n= 379, number of events= 379
##
## coef exp(coef) se(coef) z Pr(>|z|)
## airport_dist 0.125877 1.134142 0.067907 1.854 0.0638 .
## conservative -0.331611 0.717767 0.082851 -4.002 6.27e-05 ***
## male -0.141864 0.867740 0.081693 -1.737 0.0825 .
## age -0.017427 0.982724 0.103089 -0.169 0.8658
## popdens 0.108749 1.114883 0.091485 1.189 0.2346
## manufact -0.072176 0.930367 0.052325 -1.379 0.1678
## tourism 0.009962 1.010012 0.061504 0.162 0.8713
## academics -0.058010 0.943640 0.079147 -0.733 0.4636
## medinc -0.104273 0.900979 0.071349 -1.461 0.1439
## healthcare -0.173877 0.840400 0.069290 -2.509 0.0121 *
## temp 0.077350 1.080420 0.072919 1.061 0.2888
## onset_prev -0.013565 0.986527 0.075247 -0.180 0.8569
## slope_prev 0.146099 1.157311 0.086200 1.695 0.0901 .
## cpt_day_socdist_2_lag -0.238041 0.788170 0.114252 -2.083 0.0372 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## airport_dist 1.1341 0.8817 0.9928 1.2956
## conservative 0.7178 1.3932 0.6102 0.8443
## male 0.8677 1.1524 0.7394 1.0184
## age 0.9827 1.0176 0.8029 1.2028
## popdens 1.1149 0.8970 0.9319 1.3338
## manufact 0.9304 1.0748 0.8397 1.0308
## tourism 1.0100 0.9901 0.8953 1.1394
## academics 0.9436 1.0597 0.8080 1.1020
## medinc 0.9010 1.1099 0.7834 1.0362
## healthcare 0.8404 1.1899 0.7337 0.9626
## temp 1.0804 0.9256 0.9365 1.2464
## onset_prev 0.9865 1.0137 0.8513 1.1433
## slope_prev 1.1573 0.8641 0.9774 1.3703
## cpt_day_socdist_2_lag 0.7882 1.2688 0.6300 0.9860
##
## Concordance= 0.678 (se = 0.029 )
## Likelihood ratio test= 61.37 on 14 df, p=7e-08
## Wald test = 57.4 on 14 df, p=3e-07
## Score (logrank) test = 58.42 on 14 df, p=2e-07
Predict Socdist Mean
lm_mean_socdist <- spec_calculate(df = df_ger_socdist_scaled,
y = "mean_diff_socdist_2",
model = "lm",
controls = covariates %>%
append(c('mean_diff_socdist_2_lag', 'onset_prev', 'slope_prev')),
all.comb = T)
lm_mean_socdist <- lm_mean_socdist %>% map(filter_socdist)
calc_all_summaries(lm_mean_socdist)
## $pers_o
## estimate std.error statistic p.value conf.low conf.high
## mean 0.15330885 0.044355910 3.425476 0.02297009 0.06608679 0.24053090
## sd 0.07755911 0.001414087 1.656772 0.03632131 0.07585609 0.07932302
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.45206387 0.43888864 0.7476768 35.209876 7.055131e-20
## sd 0.07030338 0.07031662 0.0457296 7.949236 2.044565e-18
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 9.000000 -421.8148 865.62962 908.9425 207.11986 369.000000
## sd 1.732262 23.1611 44.31669 40.9420 26.57468 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.8122559 1 0 0.8122559
## sd 0 0.3905554 0 0 0.3905554
## significant_negative
## mean 0
## sd 0
##
## $pers_c
## estimate std.error statistic p.value conf.low conf.high
## mean -0.03851258 0.04069143 -0.9421542 0.4218812 -0.11852873 0.04150356
## sd 0.02790481 0.00308532 0.6560585 0.2885431 0.02975273 0.02730813
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.43015918 0.41651193 0.76145621 32.660794 8.760185e-07
## sd 0.09616314 0.09669053 0.06060868 9.707367 2.157079e-05
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 9.000000 -428.28432 878.56865 921.88154 215.39983 369.000000
## sd 1.732262 29.50691 57.06175 53.66655 36.34967 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.0769043 0.0546875 0.9453125 0
## sd 0 0.2664721 0.2273970 0.2273970 0
## significant_negative
## mean 0.0769043
## sd 0.2664721
##
## $pers_e
## estimate std.error statistic p.value conf.low conf.high
## mean -0.01503002 0.041886341 -0.3896659 0.5811354 -0.09739585 0.06733580
## sd 0.02278655 0.003250974 0.5237195 0.1998854 0.01851793 0.02787933
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.42867767 0.41499876 0.76240500 32.47312 3.773528e-06
## sd 0.09720623 0.09774027 0.06116183 9.72353 8.896818e-05
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 9.000000 -428.73869 879.4774 922.79028 215.95984 369.000000
## sd 1.732262 29.73164 57.4969 54.06717 36.74395 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0 0.1518555 0.8481445 0
## sd 0 0 0.3589246 0.3589246 0
## significant_negative
## mean 0
## sd 0
##
## $pers_a
## estimate std.error statistic p.value conf.low conf.high
## mean -0.02332656 0.040841977 -0.5897552 0.4122219 -0.10363874 0.05698563
## sd 0.04379788 0.002890095 1.0885441 0.2944594 0.04279858 0.04549033
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.43026173 0.41662073 0.76130558 32.712822 6.230705e-06
## sd 0.09788511 0.09844834 0.06159559 9.858796 1.500141e-04
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 9.000000 -428.17439 878.3488 921.66168 215.36107 369.000000
## sd 1.732262 29.92968 57.9119 54.51971 37.00057 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.1582031 0.3818359 0.6181641 0
## sd 0 0.3649759 0.4858960 0.4858960 0
## significant_negative
## mean 0.1582031
## sd 0.3649759
##
## $pers_n
## estimate std.error statistic p.value conf.low conf.high
## mean -0.00604797 0.041254131 -0.1365582 0.7023933 -0.08717062 0.07507468
## sd 0.01968351 0.003240101 0.4720721 0.2018550 0.02250371 0.01869854
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.42841045 0.41472473 0.76258917 32.435635 6.023785e-06
## sd 0.09718258 0.09771424 0.06110589 9.710442 1.448236e-04
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 9.000000 -428.83354 879.66708 922.97998 216.06085 369.000000
## sd 1.732262 29.68545 57.40582 53.97958 36.73502 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0 0.4108887 0.5891113 0
## sd 0 0 0.4920552 0.4920552 0
## significant_negative
## mean 0
## sd 0
plot_all_curves(lm_mean_socdist, 'lm_mean_socdist')
## $spec_results_o
##
## $spec_results_c
##
## $spec_results_e
##
## $spec_results_a
##
## $spec_results_n
lm_mean_socdist_ctrl <- lm(mean_diff_socdist_2 ~ airport_dist +
conservative + male + age + popdens + manufact +
tourism + academics + medinc + healthcare +
temp + onset_prev + slope_prev +
mean_diff_socdist_2_lag,
data = df_ger_socdist_scaled)
summary(lm_mean_socdist_ctrl)
##
## Call:
## lm(formula = mean_diff_socdist_2 ~ airport_dist + conservative +
## male + age + popdens + manufact + tourism + academics + medinc +
## healthcare + temp + onset_prev + slope_prev + mean_diff_socdist_2_lag,
## data = df_ger_socdist_scaled)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.85241 -0.44592 0.01798 0.37505 2.81498
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.001481 0.035208 -0.042 0.9665
## airport_dist -0.010038 0.045187 -0.222 0.8243
## conservative -0.273373 0.055239 -4.949 1.14e-06 ***
## male -0.117690 0.051800 -2.272 0.0237 *
## age 0.003341 0.067934 0.049 0.9608
## popdens 0.328194 0.056804 5.778 1.63e-08 ***
## manufact 0.085452 0.041025 2.083 0.0380 *
## tourism -0.013999 0.041163 -0.340 0.7340
## academics 0.364853 0.054490 6.696 8.15e-11 ***
## medinc -0.073709 0.050032 -1.473 0.1416
## healthcare -0.008339 0.046183 -0.181 0.8568
## temp -0.081000 0.048476 -1.671 0.0956 .
## onset_prev 0.008089 0.051930 0.156 0.8763
## slope_prev 0.031435 0.055945 0.562 0.5745
## mean_diff_socdist_2_lag 0.135695 0.064960 2.089 0.0374 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6853 on 364 degrees of freedom
## Multiple R-squared: 0.5478, Adjusted R-squared: 0.5304
## F-statistic: 31.49 on 14 and 364 DF, p-value: < 2.2e-16
Save all results
results_ger <- list(cox_onset_prev_hr, lm_slope_prev, lm_slope_prev_var,
lm_cases_0331, lm_cases_0930, lm_deaths_0331, lm_deaths_0930,
cox_onset_socdist_hr, lm_mean_socdist)
names(results_ger) <- list('cox_onset_prev_hr', 'lm_slope_prev', 'lm_slope_prev_var',
'lm_cases_0331', 'lm_cases_0930', 'lm_deaths_0331', 'lm_deaths_0930',
'cox_onset_socdist_hr', 'lm_mean_socdist')
save(results_ger, file="GER_spec_results.RData")
Descriptives
Distributions
df_ger_desc <- df_ger_socdist %>%
dplyr::select(pers_o, pers_c, pers_e, pers_a, pers_n,
age, male, conservative,
academics, medinc, manufact,
airport_dist, tourism, healthcare, popdens) %>%
distinct() %>% drop_na()
ger_means <- df_ger_desc %>% summarise_all(mean)
ger_sd <- df_ger_desc %>% summarise_all(sd)
desc_ger <- rbind(ger_means, ger_sd) %>% t() %>% round(3) %>% as.data.frame()
desc_ger %>% rownames_to_column() %>% write_csv('ger_descriptives.csv')
desc_ger
Explore correlations
a <- df_ger_socdist_unscaled %>% select(kreis, pers_o:pers_n, cpt_day_socdist, mean_diff_socdist)
b <- df_ger_prev_unscaled %>% select(kreis, onset_prev, slope_prev)
c <- df_ger_death_unscaled %>% select(kreis, case_rate_0331:death_rate_0930)
ger_joined <- plyr::join_all(list(a, b, c), by='kreis')
ger_joined %>% select(-kreis) %>% cor(use = 'pairwise.complete')
## pers_o pers_c pers_e pers_a pers_n
## pers_o 1.00000000 -0.10201270 0.27983871 0.15759347 -0.0640347187
## pers_c -0.10201270 1.00000000 0.21472242 0.32690021 -0.3551086811
## pers_e 0.27983871 0.21472242 1.00000000 0.18288787 -0.3649780097
## pers_a 0.15759347 0.32690021 0.18288787 1.00000000 -0.4280315714
## pers_n -0.06403472 -0.35510868 -0.36497801 -0.42803157 1.0000000000
## cpt_day_socdist -0.10416238 0.10437674 -0.04952867 -0.01672102 0.0138442293
## mean_diff_socdist 0.31429946 -0.12682312 0.09840937 -0.08341683 0.0006497013
## onset_prev -0.16574985 0.05702967 -0.27615283 -0.09470167 0.2352513229
## slope_prev 0.15470785 -0.04204164 0.24306903 0.13157797 -0.1983834025
## case_rate_0331 0.10695243 -0.01752091 0.20572870 0.13949525 -0.1722746724
## death_rate_0331 0.06566400 -0.03556727 0.10805800 0.09283184 -0.0999658180
## case_rate_0930 0.15634007 -0.05283919 0.25572247 0.04675307 -0.1116940920
## death_rate_0930 0.07587651 -0.03740828 0.08958204 0.02599047 -0.0095583322
## cpt_day_socdist mean_diff_socdist onset_prev slope_prev
## pers_o -0.1041623824 0.3142994608 -0.16574985 0.15470785
## pers_c 0.1043767434 -0.1268231184 0.05702967 -0.04204164
## pers_e -0.0495286726 0.0984093720 -0.27615283 0.24306903
## pers_a -0.0167210202 -0.0834168254 -0.09470167 0.13157797
## pers_n 0.0138442293 0.0006497013 0.23525132 -0.19838340
## cpt_day_socdist 1.0000000000 0.0358262655 0.14392128 -0.11996515
## mean_diff_socdist 0.0358262655 1.0000000000 -0.20519894 0.22948985
## onset_prev 0.1439212762 -0.2051989385 1.00000000 -0.67617865
## slope_prev -0.1199651501 0.2294898547 -0.67617865 1.00000000
## case_rate_0331 -0.0689721829 0.1733238658 -0.66024501 0.92767374
## death_rate_0331 0.0003684063 0.1435286041 -0.40795429 0.69587384
## case_rate_0930 -0.1260015001 0.2185131719 -0.58477976 0.82060844
## death_rate_0930 -0.0078242968 0.1456354879 -0.32697005 0.64249969
## case_rate_0331 death_rate_0331 case_rate_0930 death_rate_0930
## pers_o 0.10695243 0.0656640020 0.15634007 0.075876508
## pers_c -0.01752091 -0.0355672682 -0.05283919 -0.037408282
## pers_e 0.20572870 0.1080579960 0.25572247 0.089582044
## pers_a 0.13949525 0.0928318364 0.04675307 0.025990472
## pers_n -0.17227467 -0.0999658180 -0.11169409 -0.009558332
## cpt_day_socdist -0.06897218 0.0003684063 -0.12600150 -0.007824297
## mean_diff_socdist 0.17332387 0.1435286041 0.21851317 0.145635488
## onset_prev -0.66024501 -0.4079542869 -0.58477976 -0.326970047
## slope_prev 0.92767374 0.6958738421 0.82060844 0.642499692
## case_rate_0331 1.00000000 0.7967609328 0.82553862 0.687015730
## death_rate_0331 0.79676093 1.0000000000 0.68600498 0.850162063
## case_rate_0930 0.82553862 0.6860049833 1.00000000 0.745023167
## death_rate_0930 0.68701573 0.8501620632 0.74502317 1.000000000
ger_means <- ger_joined %>% select(-kreis) %>% summarise_all(mean)
ger_sd <- ger_joined %>% select(-kreis) %>% summarise_all(sd)
desc_ger <- rbind(ger_means, ger_sd) %>% t() %>% round(3) %>% as.data.frame() %>%
rename(mean=V1, sd=V2)
desc_ger %>% rownames_to_column() %>% write_csv('ger_descriptives.csv')
desc_ger
Visualizations
df_ger_prev <- read_csv('GER_prevalence.csv')
df_ger_prev <- df_ger_prev %>%
mutate(date = as.Date(date, "%d%b%Y"),
kreis = as.character(kreis)) %>%
dplyr::select(kreis, date, rate_day) %>%
filter(date < '2020-04-01')
# create sequence of dates
date_sequence <- seq.Date(min(df_ger_prev$date),
max(df_ger_prev$date), 1)
# create data frame with time sequence
df_dates = data.frame(date_sequence, 1:length(date_sequence))
names(df_dates) <- c('date', 'time')
tto <- df_ger_prev_unscaled %>%
filter(event==1) %>%
merge(df_dates, by.x = 'onset_prev', by.y = 'time') %>%
ggplot(aes(x=date+1)) +
geom_histogram(bins=15) +
theme_light() +
ggtitle('Distribution - COVID-19 Onset') +
theme(plot.title = element_text(size=11)) +
ylab('Absolute frequency') +
xlab('Date')
soa <- df_ger_slope_prev %>%
ggplot(aes(x=slope_prev)) +
geom_histogram(bins=15) +
theme_light() +
ggtitle('Distribution - COVID-19 Growth Rates') +
theme(plot.title = element_text(size=11)) +
ylab('Absolute frequency') +
xlab('Standardized growth rates')
figure <- ggarrange(tto, soa,
labels = c("", ""),
ncol = 2, nrow = 1)
figure
tto <- df_ger_death %>%
ggplot(aes(x=case_rate_0930)) +
geom_histogram(bins=15) +
theme_light() +
ggtitle('Distribution - Cumulative Case Rates') +
theme(plot.title = element_text(size=11)) +
ylab('Absolute frequency') +
xlab('Cases Per 1000 Inhabitants')
soa <- df_ger_death %>%
ggplot(aes(x=death_rate_0930)) +
geom_histogram(bins=15) +
theme_light() +
ggtitle('Distribution - Cumulative Death Rates') +
theme(plot.title = element_text(size=11)) +
ylab('Absolute frequency') +
xlab('Deaths Per 1000 Inhabitants')
figure <- ggarrange(tto, soa,
labels = c("", ""),
ncol = 2, nrow = 1)
figure
# create sequence of dates
date_sequence <- seq.Date(min(df_ger_socdist$date),
max(df_ger_socdist$date), 1)
# create data frame with time sequence
df_dates = data.frame(date_sequence, 1:length(date_sequence))
names(df_dates) <- c('date', 'time')
tto <- df_ger_cpt_socdist %>%
merge(df_dates, by.x = 'cpt_day_socdist_2', by.y = 'time') %>%
#filter(cpt_day_socdist_2 < 30) %>%
ggplot(aes(x=date+1)) +
geom_histogram(bins=15) +
theme_light() +
ggtitle('Distribution - Time to Adoption') +
theme(plot.title = element_text(size=11)) +
ylab('Absolute frequency') +
xlab('Date') +
xlim(c(as.Date('2020-03-01'), as.Date('2020-03-31')))
soa <- df_ger_cpt_socdist %>%
ggplot(aes(x=mean_diff_socdist_2)) +
geom_histogram(bins=15) +
theme_light() +
ggtitle('Distribution - Strength of Adjustment') +
theme(plot.title = element_text(size=11)) +
ylab('Absolute frequency') +
xlab('Standardized mean difference')
figure <- ggarrange(tto, soa,
labels = c("", ""),
ncol = 2, nrow = 1)
figure
=======
Prepare data
Merge prevalence data
# create sequence of dates
date_sequence <- seq.Date(min(df_ger_prev$date),
max(df_ger_prev$date), 1)
# create data frame with time sequence
df_dates = tibble(date_sequence, 1:length(date_sequence))
names(df_dates) <- c('date', 'time')
# merge prevalence data
df_ger_prev <- df_ger_prev %>%
inner_join(df_ger_pers, by = 'kreis') %>%
inner_join(df_ger_ctrl, by = 'kreis') %>%
merge(df_dates, by='date') %>%
arrange(kreis, date) %>%
drop_na()
df_ger_prev %>% select(kreis) %>% distinct() %>% nrow()
## [1] 383
# join data frames
df_ger_death <- df_ger_death %>%
plyr::join(df_ger_ctrl, by='kreis') %>%
plyr::join(df_ger_pers, by='kreis') %>%
drop_na()
# save df
df_ger_death %>%
#select(kreis, case_rate, death_rate) %>%
write_csv('ger_death_kreis_controls.csv')
Merge social distancing data
# create sequence of dates
date_sequence <- seq.Date(min(df_ger_socdist$date),
max(df_ger_socdist$date), 1)
# create data frame with time sequence
df_dates = tibble(date_sequence, 1:length(date_sequence))
names(df_dates) <- c('date', 'time')
# merge socdist data
df_ger_socdist <- df_ger_socdist %>%
inner_join(df_ger_pers, by = 'kreis') %>%
inner_join(df_ger_ctrl, by = 'kreis') %>%
merge(df_dates, by='date') %>%
arrange(kreis)
df_ger_socdist %>% select(kreis) %>% distinct() %>% nrow()
## [1] 379
Control for weekend effect
easter <- seq.Date(as.Date('2020-04-10'), as.Date('2020-04-13'), 1)
df_ger_loess <- df_ger_socdist %>%
mutate(weekday = format(date, '%u')) %>%
filter(!(weekday %in% c('6','7') | date %in% easter)) %>%
split(.$kreis) %>%
map(~ loess(socdist_single_tile ~ time, data = .)) %>%
map(predict, 1:max(df_ger_socdist$time)) %>%
bind_rows() %>%
gather(key = 'kreis', value = 'loess') %>%
group_by(kreis) %>%
mutate(time = row_number())
df_ger_loess_2 <- df_ger_socdist %>%
mutate(weekday = format(date, '%u')) %>%
filter(!(weekday %in% c('6','7') | date %in% easter)) %>%
split(.$kreis) %>%
map(~ loess(socdist_tiles ~ time, data = .)) %>%
map(predict, 1:max(df_ger_socdist$time)) %>%
bind_rows() %>%
gather(key = 'kreis', value = 'loess') %>%
rename(loess_2 = loess) %>%
group_by(kreis) %>%
mutate(time = row_number())
df_ger_socdist <- df_ger_socdist %>%
merge(df_ger_loess, by=c('kreis', 'time')) %>%
merge(df_ger_loess_2, by=c('kreis', 'time')) %>%
mutate(weekday = format(date, '%u')) %>%
mutate(socdist_single_tile_clean = ifelse(weekday %in% c('6','7') | date %in% easter,
loess, socdist_single_tile),
socdist_tiles_clean = ifelse(weekday %in% c('6','7') | date %in% easter,
loess_2, socdist_tiles)) %>%
arrange(kreis, time) %>%
select(-weekday)
df_ger_socdist <- df_ger_socdist %>% drop_na() %>% mutate(time = time-1)
Include cleaned social distancing data
df_ger_socdist <- df_ger_socdist %>%
mutate(socdist_single_tile = socdist_single_tile_clean,
socdist_tiles = socdist_tiles_clean) %>%
select(-loess, -loess_2, -socdist_single_tile_clean, -socdist_tiles_clean)
Define outcomes of interest
Social Distancing - Change point analysis
# keep only counties with full data
kreis_complete <- df_ger_socdist %>%
group_by(kreis) %>%
summarise(n = n()) %>%
filter(n==max(.$n)) %>%
.$kreis
# run changepoint analysis
df_ger_socdist_cpt_results <- df_ger_socdist %>%
select(kreis, socdist_single_tile) %>%
filter(kreis %in% kreis_complete) %>%
split(.$kreis) %>%
map(~ cpt.meanvar(as.vector(.$socdist_single_tile),
#penalty = 'Asymptotic',
class=TRUE,
param.estimates=TRUE,
Q=1,
test.stat = 'Normal'))
df_ger_socdist_cpt_results_2 <- df_ger_socdist %>%
select(kreis, socdist_tiles) %>%
filter(kreis %in% kreis_complete) %>%
split(.$kreis) %>%
map(~ cpt.meanvar(as.vector(.$socdist_tiles),
class=TRUE,
param.estimates=TRUE,
Q=1,
test.stat = 'Normal'))
# calculate change point
df_ger_socdist_cpt_day <- df_ger_socdist_cpt_results %>%
map(cpts) %>%
unlist() %>%
as.data.frame() %>%
rename(cpt_day_socdist = '.') %>%
rownames_to_column('kreis')
df_ger_socdist_cpt_day_2 <- df_ger_socdist_cpt_results_2 %>%
map(cpts) %>%
unlist() %>%
as.data.frame() %>%
rename(cpt_day_socdist_2 = '.') %>%
rownames_to_column('kreis')
# calculate mean differences
df_ger_socdist_cpt_mean_diff <- df_ger_socdist_cpt_results %>%
map(param.est) %>%
map(~ .$mean) %>%
map(~ .[2]) %>%
unlist() %>%
as.data.frame() %>%
rename(mean_diff_socdist = '.') %>%
rownames_to_column('kreis')
df_ger_socdist_cpt_mean_diff_2 <- df_ger_socdist_cpt_results_2 %>%
map(param.est) %>%
map(~ .$mean) %>%
map(~ -.[2]) %>%
unlist() %>%
as.data.frame() %>%
rename(mean_diff_socdist_2 = '.') %>%
rownames_to_column('kreis')
# calculate means
df_ger_socdist_mean <- df_ger_socdist %>%
group_by(kreis) %>%
summarise(mean_socdist = mean(socdist_single_tile))
df_ger_socdist_mean_2 <- df_ger_socdist %>%
group_by(kreis) %>%
summarise(mean_socdist_2 = -mean(socdist_tiles))
# merge with county data
df_ger_cpt_socdist <- df_ger_socdist %>%
select(-date, -time, -socdist_single_tile, -socdist_tiles) %>%
distinct() %>%
left_join(df_ger_socdist_cpt_day, by='kreis') %>%
left_join(df_ger_socdist_cpt_day_2, by='kreis') %>%
left_join(df_ger_socdist_cpt_mean_diff, by='kreis') %>%
left_join(df_ger_socdist_cpt_mean_diff_2, by='kreis') %>%
left_join(df_ger_socdist_mean, by='kreis') %>%
left_join(df_ger_socdist_mean_2, by='kreis') %>%
left_join(select(df_ger_onset_prev, kreis, onset_prev), by='kreis') %>%
left_join(select(df_ger_slope_prev, kreis, slope_prev), by='kreis')
# handle censored data
df_ger_cpt_socdist <- df_ger_cpt_socdist %>%
mutate(cpt_day_socdist = ifelse(is.na(cpt_day_socdist),
as.numeric(diff(range(df_ger_socdist$date))),
cpt_day_socdist)) %>%
mutate(event = ifelse(cpt_day_socdist >=
as.numeric(diff(range(df_ger_socdist$date))), 0, 1))
Remove incomplete cases
df_ger_prev_unscaled <- inner_join(df_ger_slope_prev, df_ger_cpt_socdist %>% select(kreis), by='kreis')
df_ger_death_unscaled <- inner_join(df_ger_death, df_ger_cpt_socdist %>% select(kreis), by='kreis')
df_ger_socdist_unscaled <- inner_join(df_ger_cpt_socdist, df_ger_slope_prev %>% select(kreis), by='kreis')
kreis_list <- df_ger_prev_unscaled %>% select(kreis)
kreis_list %>% write_csv('kreis_list.csv')
df_ger_prev_unscaled %>% merge(df_ger_death_unscaled %>% select(kreis:death_rate_0930)) %>%
merge(df_ger_socdist_unscaled %>% select(kreis, cpt_day_socdist_2, mean_diff_socdist_2)) %>%
select(-event) %>% write_csv('ger_all_variables.csv')
Rescale Data
df_ger_prev_scaled <- df_ger_prev_unscaled %>%
mutate_at(vars(-kreis, -event), scale)
df_ger_death_scaled <- df_ger_death_unscaled %>%
mutate_at(vars(-kreis), scale)
df_ger_socdist_scaled <- df_ger_socdist_unscaled %>%
mutate_at(vars(-kreis, -event), scale)
Calculate spatial lags
# create UDF to calculate lagged variables
calc_lags <- function(df, weights, cols){
cols_only <- df %>% select(cols)
cols_lag <- weights %*% as.matrix(cols_only) %>%
as.matrix() %>% as.data.frame()
names(cols_lag) <- paste0(names(cols_lag), '_lag')
return(cols_lag)
}
# read_weight matrix
weight_mat_norm <- read_csv('GER_spatial_weights.csv')
weight_mat_norm <- weight_mat_norm %>% select(-kreis) %>% as.matrix()
dim(weight_mat_norm)
## [1] 379 379
# generate spatially lagged y
y_only <- df_ger_prev_scaled %>% select(onset_prev,slope_prev,slope_prev_var) %>% names()
y_lag <- calc_lags(df_ger_prev_scaled, weight_mat_norm, y_only)
# generate spatially lagged X
X_only <- df_ger_prev_scaled %>% select(academics:pers_n) %>% names()
X_lag <- calc_lags(df_ger_prev_scaled, weight_mat_norm, X_only)
# bind new variables to df
df_ger_prev_scaled <- cbind(df_ger_prev_scaled, y_lag, X_lag)
# generate spatially lagged y
y_only <- df_ger_death_scaled %>% select(case_rate_0331:death_rate_0930) %>% names()
y_lag <- calc_lags(df_ger_death_scaled, weight_mat_norm, y_only)
# generate spatially lagged X
X_only <- df_ger_death_scaled %>% select(academics:pers_n) %>% names()
X_lag <- calc_lags(df_ger_death_scaled, weight_mat_norm, X_only)
# bind new variables to df
df_ger_death_scaled <- cbind(df_ger_death_scaled, y_lag, X_lag)
# generate spatially lagged y
y_only <- df_ger_socdist_scaled %>% select(cpt_day_socdist_2,mean_diff_socdist_2) %>% names()
y_lag <- calc_lags(df_ger_socdist_scaled, weight_mat_norm, y_only)
# generate spatially lagged X
X_only <- df_ger_socdist_scaled %>% select(academics:pers_n) %>% names()
X_lag <- calc_lags(df_ger_socdist_scaled, weight_mat_norm, X_only)
# bind new variables to df
df_ger_socdist_scaled <- cbind(df_ger_socdist_scaled, y_lag, X_lag)
write_csv(df_ger_prev_scaled, 'df_ger_slope_prev.csv')
write_csv(df_ger_death_scaled, 'df_ger_death_scaled.csv')
write_csv(df_ger_socdist_scaled, 'df_ger_cpt_socdist.csv')
Run Specification Curve Analysis
# define function to calculate specification curve analyses for all traits
spec_calculate <- function(df, y, model, controls, all.comb = T){
spec_results_o <- run_specs(df = df,
y = c(y), x = c("pers_o"),
model = c(model), controls = controls,
all.comb = all.comb)
spec_results_c <- run_specs(df = df,
y = c(y), x = c("pers_c"),
model = c(model), controls = controls,
all.comb = all.comb)
spec_results_e <- run_specs(df = df,
y = c(y), x = c("pers_e"),
model = c(model), controls = controls,
all.comb = all.comb)
spec_results_a <- run_specs(df = df,
y = c(y), x = c("pers_a"),
model = c(model), controls = controls,
all.comb = all.comb)
spec_results_n <- run_specs(df = df,
y = c(y), x = c("pers_n"),
model = c(model), controls = controls,
all.comb = all.comb)
spec_results <- list(spec_results_o, spec_results_c, spec_results_e,
spec_results_a, spec_results_n)
names(spec_results) <- list('spec_results_o', 'spec_results_c', 'spec_results_e',
'spec_results_a', 'spec_results_n')
return(spec_results)
}
# adapted from specr package code
format_results <- function(df, var, null = 0, desc = FALSE) {
# rank specs
if (isFALSE(desc)) {
df <- df %>%
dplyr::arrange(!! var)
} else {
df <- df %>%
dplyr::arrange(desc(!! var))
}
# create rank variable and color significance
df <- df %>%
dplyr::mutate(specifications = 1:nrow(df),
color = case_when(conf.low > null ~ "#377eb8",
conf.high < null ~ "#e41a1c",
TRUE ~ "darkgrey"))
return(df)
}
# define function to plot single specification curve
plot_curves <- function(results, filename, hr=F){
file_path_eps <- paste0("../../Plots/GER/", filename,".eps")
file_path_pdf <- paste0("../../Plots/GER/", filename,".pdf")
file_path_png <- paste0("../../Plots/GER/", filename,".png")
if(hr==F){
results %>%
format_results(var = .$estimate, null = 0, desc = F) %>%
ggplot(aes(x = specifications,
y = estimate,
ymin = conf.low,
ymax = conf.high,
color = color)) +
geom_point(aes(color = color),
size = 1) +
theme_minimal() +
scale_color_identity() +
theme(strip.text = element_blank(),
axis.line = element_line("black", size = .5),
legend.position = "none",
panel.spacing = unit(.75, "lines"),
axis.text = element_text(colour = "black")) +
labs(x = "") +
geom_pointrange(alpha = 0.05,
size = .6,
fatten = 1) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
labs(x = "", y = "standarized coefficient") +
coord_fixed(ratio = 2000, ylim = c(-0.4, 0.4)) +
theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black")) +
ggsave(file=file_path_eps, device = 'eps')+
ggsave(file=file_path_pdf, device = 'pdf')+
ggsave(file=file_path_png, device = 'png')
}else{
results %>%
format_results(var = .$estimate, null = 1, desc = F) %>%
ggplot(aes(x = specifications,
y = estimate,
ymin = conf.low,
ymax = conf.high,
color = color)) +
geom_point(aes(color = color),
size = 1) +
theme_minimal() +
scale_color_identity() +
theme(strip.text = element_blank(),
axis.line = element_line("black", size = .5),
legend.position = "none",
panel.spacing = unit(.75, "lines"),
axis.text = element_text(colour = "black")) +
labs(x = "") +
geom_pointrange(alpha = 0.05,
size = .6,
fatten = 1) +
geom_hline(yintercept = 1, linetype = "dashed", color = "black") +
labs(x = "", y = "standarized coefficient") +
coord_fixed(ratio = 2000, ylim = c(0.6, 1.4)) +
theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black")) +
ggsave(file=file_path_eps, device = 'eps')+
ggsave(file=file_path_pdf, device = 'pdf')+
ggsave(file=file_path_png, device = 'png')
}
}
# define function to plot all curves in list
plot_all_curves <- function(ls, analysis, hr=F){
pers <- c('o', 'c', 'e', 'a', 'n')
filenames <- as.list(paste0(analysis, '_', pers))
hr <- as.list(rep(hr, 5))
pmap(list(ls, filenames, hr), plot_curves)
}
# define function to calculate summary stats for single specification curve
calc_summary <- function(df){
dft <- df %>% select(estimate:fit_nobs)
dft <- dft %>% mutate(significant = as.numeric(p.value<0.05),
positive = as.numeric(statistic>0),
negative = as.numeric(statistic<0),
significant_positive = as.numeric(p.value<0.05 & statistic>0),
significant_negative = as.numeric(p.value<0.05 & statistic<0))
mean_temp <- dft %>% map_if(is.numeric, mean, .else=NULL) %>% as.data.frame()
sd_temp <- dft %>% map_if(is.numeric, sd, .else=NULL) %>% as.data.frame()
df_temp <- rbind(mean_temp, sd_temp)
row.names(df_temp) <- c('mean', 'sd')
return(df_temp)
}
# define function to calculate summary stats for all curves in list
calc_all_summaries <- function(ls){
ls <- ls %>% map(calc_summary)
names(ls) <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')
return(ls)
}
coef_to_hr <- function(df){
df %>% mutate(conf.low = exp(estimate-1.96*std.error),
conf.high = exp(estimate+1.96*std.error),
estimate = exp(estimate))
}
filter_socdist <- function(df){
df %>% filter(str_detect(controls, 'onset_prev') &
str_detect(controls, 'slope_prev'))
}
covariates <- c("airport_dist", "conservative", 'male', 'age', 'popdens',
'manufact', 'tourism', 'academics', 'medinc', 'healthcare',
'temp')
Predict Onset
cox_model <- function(formula, data){
formula <- as.formula(formula)
coxph(formula = formula, data = data)}
cox_onset_prev <- spec_calculate(df = df_ger_prev_scaled,
y = "Surv(onset_prev, event)",
model = "cox_model",
controls = covariates %>% append('onset_prev_lag'),
all.comb = T)
cox_onset_prev_hr <- cox_onset_prev %>% map(coef_to_hr)
calc_all_summaries(cox_onset_prev_hr)
## $pers_o
## estimate std.error statistic p.value conf.low conf.high fit_n
## mean 1.02726691 0.055320547 0.4779412 0.5231112 0.92180104 1.14483537 379
## sd 0.05390301 0.002840416 0.9535330 0.2800335 0.05018495 0.05845821 0
## fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean 379 131.39127 3.811112e-06 140.8145
## sd 0 36.26568 5.698238e-05 41.9226
## fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean 4.156490e-06 137.52997 4.234893e-06 NA
## sd 5.770949e-05 40.56007 5.819916e-05 NA
## fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean NA 0.2895930 0.9999496 0.7195157
## sd NA 0.0722395 0.0000000 0.0381976
## fit_std.error.concordance fit_logLik fit_AIC fit_BIC fit_nobs
## mean 0.015635755 -1809.51851 3633.0370 3660.59978 379
## sd 0.001317056 18.13284 34.3554 31.38233 0
## significant positive negative significant_positive significant_negative
## mean 0.08911133 0.6496582 0.3503418 0.08911133 0
## sd 0.28493915 0.4771352 0.4771352 0.28493915 0
##
## $pers_c
## estimate std.error statistic p.value conf.low conf.high fit_n
## mean 0.94922252 0.053595445 -0.9954995 0.3479716 0.85456571 1.05437271 379
## sd 0.04503002 0.001467328 0.8933022 0.2887067 0.04054931 0.05018395 0
## fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean 379 132.01667 0.0001954177 140.8434
## sd 0 37.85635 0.0029957287 43.1195
## fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean 0.0001852582 137.59893 0.0001881611 NA
## sd 0.0028016025 41.78484 0.0028121428 NA
## fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean NA 0.29044877 0.9999496 0.71826746
## sd NA 0.07571748 0.0000000 0.04121107
## fit_std.error.concordance fit_logLik fit_AIC fit_BIC fit_nobs
## mean 0.015744358 -1809.20581 3632.41162 3659.9744 379
## sd 0.001364156 18.92817 35.97484 33.0281 0
## significant positive negative significant_positive significant_negative
## mean 0.1445312 0.1584473 0.8415527 0 0.1445312
## sd 0.3516705 0.3652045 0.3652045 0 0.3516705
##
## $pers_e
## estimate std.error statistic p.value conf.low conf.high fit_n
## mean 1.10484855 0.052014824 1.9131109 0.1418816 0.99784541 1.22333370 379
## sd 0.05117709 0.001260917 0.9230424 0.1475011 0.04805772 0.05449954 0
## fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean 379 134.80506 6.667657e-08 144.44569
## sd 0 33.69771 7.667802e-07 39.58589
## fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean 1.086854e-07 141.00631 1.184551e-07 NA
## sd 1.098564e-06 38.16651 1.180765e-06 NA
## fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean NA 0.29642439 0.9999496 0.72680567
## sd NA 0.06635634 0.0000000 0.03214107
## fit_std.error.concordance fit_logLik fit_AIC fit_BIC fit_nobs
## mean 0.015317233 -1807.81162 3629.62324 3657.18599 379
## sd 0.001142571 16.84885 31.78596 28.87968 0
## significant positive negative significant_positive significant_negative
## mean 0.3879395 1 0 0.3879395 0
## sd 0.4873401 0 0 0.4873401 0
##
## $pers_a
## estimate std.error statistic p.value conf.low conf.high fit_n
## mean 1.006651 0.053837899 0.1292402 0.645567 0.90592390 1.11859211 379
## sd 0.032338 0.001749636 0.6071196 0.238925 0.03165915 0.03294262 0
## fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean 379 130.6175 0.0006586847 140.37322
## sd 0 37.6386 0.0123791555 43.20648
## fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean 0.0006630851 137.15705 0.000668293 NA
## sd 0.0122500086 41.79493 0.012237339 NA
## fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean NA 0.2878630 0.9999496 0.71879483
## sd NA 0.0756924 0.0000000 0.04050933
## fit_std.error.concordance fit_logLik fit_AIC fit_BIC fit_nobs
## mean 0.015616239 -1809.9054 3633.8108 3661.37351 379
## sd 0.001276486 18.8193 35.7401 32.76038 0
## significant positive negative significant_positive significant_negative
## mean 0 0.4870605 0.5129395 0 0
## sd 0 0.4998936 0.4998936 0 0
##
## $pers_n
## estimate std.error statistic p.value conf.low conf.high fit_n
## mean 0.8618104 0.0536108475 -2.7813176 0.01334193 0.77585634 0.95728861 379
## sd 0.0256510 0.0006767025 0.5515346 0.01683804 0.02333785 0.02824371 0
## fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean 379 138.42138 8.297546e-07 148.34185
## sd 0 35.84265 1.236878e-05 41.64845
## fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean 9.623814e-07 144.57863 1.004892e-06 NA
## sd 1.340854e-05 40.36034 1.343165e-05 NA
## fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean NA 0.30271384 0.9999496 0.72815172
## sd NA 0.07048457 0.0000000 0.03346175
## fit_std.error.concordance fit_logLik fit_AIC fit_BIC fit_nobs
## mean 0.015109031 -1806.00346 3626.00692 3653.56967 379
## sd 0.001124059 17.92132 33.95795 31.05319 0
## significant positive negative significant_positive significant_negative
## mean 0.9492188 0 1 0 0.9492188
## sd 0.2195775 0 0 0 0.2195775
plot_all_curves(cox_onset_prev_hr, 'cox_onset_prev_hr', hr = T)
## $spec_results_o
##
## $spec_results_c
##
## $spec_results_e
##
## $spec_results_a
##
## $spec_results_n
cox_onset_prev_ctrl <- coxph(Surv(onset_prev, event) ~ airport_dist +
conservative + male + age + popdens + manufact +
tourism + academics + medinc + healthcare +
temp + onset_prev_lag,
data = df_ger_prev_scaled)
summary(cox_onset_prev_ctrl)
## Call:
## coxph(formula = Surv(onset_prev, event) ~ airport_dist + conservative +
## male + age + popdens + manufact + tourism + academics + medinc +
## healthcare + temp + onset_prev_lag, data = df_ger_prev_scaled)
##
## n= 379, number of events= 379
##
## coef exp(coef) se(coef) z Pr(>|z|)
## airport_dist -0.09147 0.91259 0.07045 -1.298 0.194133
## conservative -0.08051 0.92264 0.08006 -1.006 0.314608
## male -0.13923 0.87003 0.07451 -1.869 0.061664 .
## age -0.25167 0.77750 0.09118 -2.760 0.005779 **
## popdens -0.09156 0.91251 0.07416 -1.235 0.216965
## manufact 0.00861 1.00865 0.05321 0.162 0.871458
## tourism -0.11689 0.88968 0.05988 -1.952 0.050933 .
## academics 0.07279 1.07551 0.08113 0.897 0.369583
## medinc 0.29682 1.34558 0.07898 3.758 0.000171 ***
## healthcare -0.06525 0.93683 0.06938 -0.940 0.347005
## temp -0.07546 0.92731 0.07139 -1.057 0.290488
## onset_prev_lag -0.36945 0.69111 0.08959 -4.124 3.73e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## airport_dist 0.9126 1.0958 0.7949 1.0477
## conservative 0.9226 1.0838 0.7886 1.0794
## male 0.8700 1.1494 0.7518 1.0068
## age 0.7775 1.2862 0.6503 0.9296
## popdens 0.9125 1.0959 0.7891 1.0553
## manufact 1.0086 0.9914 0.9087 1.1195
## tourism 0.8897 1.1240 0.7912 1.0005
## academics 1.0755 0.9298 0.9174 1.2609
## medinc 1.3456 0.7432 1.1526 1.5709
## healthcare 0.9368 1.0674 0.8177 1.0733
## temp 0.9273 1.0784 0.8062 1.0666
## onset_prev_lag 0.6911 1.4469 0.5798 0.8238
##
## Concordance= 0.755 (se = 0.015 )
## Likelihood ratio test= 177.8 on 12 df, p=<2e-16
## Wald test = 183.9 on 12 df, p=<2e-16
## Score (logrank) test = 193.2 on 12 df, p=<2e-16
Predict Slopes
# fixed time windows
lm_slope_prev <- spec_calculate(df = df_ger_prev_scaled,
y = "slope_prev",
model = "lm",
controls = covariates %>% append('slope_prev_lag'),
all.comb = T)
calc_all_summaries(lm_slope_prev)
## $pers_o
## estimate std.error statistic p.value conf.low conf.high
## mean 0.04770701 0.047105679 0.9732636 0.4496175 -0.04492057 0.14033459
## sd 0.04866155 0.003961232 0.9445598 0.3249956 0.04463451 0.05352556
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3762331 0.3647529 0.79419554 34.39651 0.0000103227
## sd 0.1110695 0.1112149 0.06709402 12.97273 0.0001987136
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -445.11260 908.22521 943.66303 235.78389 371.000000
## sd 1.732262 31.40818 60.85522 57.41092 41.98425 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.1411133 0.8693848 0.1306152 0.1411133
## sd 0 0.3481809 0.3370202 0.3370202 0.3481809
## significant_negative
## mean 0
## sd 0
##
## $pers_c
## estimate std.error statistic p.value conf.low conf.high
## mean -0.02472864 0.042382347 -0.5866914 0.4324615 -0.10806835 0.05861107
## sd 0.04143443 0.003375866 0.9654986 0.2623563 0.04195917 0.04196624
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3747063 0.3632077 0.79495961 34.26618 0.001065338
## sd 0.1157000 0.1159113 0.06951753 13.36336 0.018659191
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -445.39270 908.78539 944.22322 236.36101 371.000000
## sd 1.732262 32.35788 62.76145 59.31843 43.73461 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.09521484 0.3103027 0.6896973 0
## sd 0 0.29354730 0.4626740 0.4626740 0
## significant_negative
## mean 0.09521484
## sd 0.29354730
##
## $pers_e
## estimate std.error statistic p.value conf.low conf.high
## mean 0.09475842 0.042916475 2.145207 0.1152644 0.01036841 0.17914843
## sd 0.05415151 0.002851794 1.072831 0.1384433 0.04912212 0.05928462
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3835310 0.3721721 0.78983225 35.33176 6.972021e-08
## sd 0.1038293 0.1038635 0.06319739 12.45035 1.032508e-06
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -443.15352 904.30703 939.74486 233.02529 371.000000
## sd 1.732262 29.85157 57.72422 54.26548 39.24747 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.4655762 1 0 0.4655762
## sd 0 0.4988745 0 0 0.4988745
## significant_negative
## mean 0
## sd 0
##
## $pers_a
## estimate std.error statistic p.value conf.low conf.high
## mean 0.04444294 0.042531589 1.0138143 0.4314829 -0.03919023 0.12807612
## sd 0.03737716 0.003433074 0.8010556 0.3214330 0.03396617 0.04161167
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3757420 0.3642556 0.7944327 34.36378 0.0001083107
## sd 0.1128081 0.1129825 0.0679861 13.13045 0.0019291217
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -445.19538 908.3908 943.82859 235.96953 371.000000
## sd 1.732262 31.75005 61.5464 58.11313 42.64145 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.1437988 0.9528809 0.04711914 0.1437988
## sd 0 0.3509285 0.2119195 0.21191954 0.3509285
## significant_negative
## mean 0
## sd 0
##
## $pers_n
## estimate std.error statistic p.value conf.low conf.high
## mean -0.06378630 0.042452946 -1.458345 0.3027720 -0.14726483 0.01969223
## sd 0.04686406 0.003144902 0.979003 0.2759218 0.05132561 0.04283259
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3785076 0.3670643 0.79281621 34.71477 2.890617e-06
## sd 0.1093210 0.1094544 0.06617581 12.90910 5.052404e-05
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -444.48386 906.96772 942.40555 234.92414 371.000000
## sd 1.732262 31.04716 60.14311 56.72546 41.32336 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.3442383 0 1 0
## sd 0 0.4751772 0 0 0
## significant_negative
## mean 0.3442383
## sd 0.4751772
plot_all_curves(lm_slope_prev, 'lm_slope_prev')
## $spec_results_o
##
## $spec_results_c
##
## $spec_results_e
##
## $spec_results_a
##
## $spec_results_n
lm_slope_prev_ctrl <- lm(slope_prev ~ airport_dist +
conservative + male + age + popdens + manufact +
tourism + academics + medinc + healthcare +
temp + slope_prev_lag,
data = df_ger_prev_scaled)
summary(lm_slope_prev_ctrl)
##
## Call:
## lm(formula = slope_prev ~ airport_dist + conservative + male +
## age + popdens + manufact + tourism + academics + medinc +
## healthcare + temp + slope_prev_lag, data = df_ger_prev_scaled)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.5022 -0.4555 -0.1231 0.3315 3.8222
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.007574 0.036347 -0.208 0.835056
## airport_dist -0.010112 0.046450 -0.218 0.827784
## conservative 0.007307 0.056396 0.130 0.896981
## male -0.085540 0.053273 -1.606 0.109207
## age -0.296574 0.069586 -4.262 2.58e-05 ***
## popdens -0.116212 0.058141 -1.999 0.046368 *
## manufact 0.061335 0.042127 1.456 0.146268
## tourism -0.126881 0.041916 -3.027 0.002645 **
## academics 0.025668 0.056486 0.454 0.649799
## medinc 0.239428 0.052239 4.583 6.29e-06 ***
## healthcare -0.018543 0.047757 -0.388 0.698034
## temp -0.174761 0.049146 -3.556 0.000426 ***
## slope_prev_lag 0.432857 0.060688 7.133 5.31e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7073 on 366 degrees of freedom
## Multiple R-squared: 0.5156, Adjusted R-squared: 0.4997
## F-statistic: 32.47 on 12 and 366 DF, p-value: < 2.2e-16
# variable time windows
lm_slope_prev_var <- spec_calculate(df = df_ger_prev_scaled,
y = "slope_prev_var",
model = "lm",
controls = covariates %>% append('slope_prev_var_lag'),
all.comb = T)
calc_all_summaries(lm_slope_prev_var)
## $pers_o
## estimate std.error statistic p.value conf.low conf.high
## mean 0.03138877 0.048027782 0.6168286 0.5240564 -0.06305202 0.1258296
## sd 0.04551793 0.003904359 0.8881225 0.3037374 0.04165726 0.0502624
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3513107 0.3393623 0.81020283 30.71510 0.0001379822
## sd 0.1089591 0.1089781 0.06488532 11.75029 0.0027232807
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -452.79671 923.59341 959.03124 245.20455 371.000000
## sd 1.732262 29.99134 57.97282 54.44501 41.18655 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.07739258 0.7424316 0.2575684 0.07739258
## sd 0 0.26724596 0.4373484 0.4373484 0.26724596
## significant_negative
## mean 0
## sd 0
##
## $pers_c
## estimate std.error statistic p.value conf.low conf.high
## mean -0.005412083 0.043251658 -0.1122928 0.4938687 -0.09046118 0.07963702
## sd 0.041568173 0.003203958 0.9510700 0.2765537 0.04319003 0.04086338
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3505259 0.3385690 0.81057266 30.65762 0.001475911
## sd 0.1118288 0.1118893 0.06636315 11.99597 0.025347073
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -452.92184 923.8437 959.28150 245.50122 371.000000
## sd 1.732262 30.55817 59.1121 55.58766 42.27129 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.04272461 0.4780273 0.5219727 0
## sd 0 0.20226024 0.4995780 0.4995780 0
## significant_negative
## mean 0.04272461
## sd 0.20226024
##
## $pers_e
## estimate std.error statistic p.value conf.low conf.high
## mean 0.06672988 0.043867749 1.464399 0.3014781 -0.01953070 0.1529904
## sd 0.05252357 0.002821869 1.068661 0.2892869 0.04787854 0.0573294
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3555855 0.3437053 0.80774108 31.20591 2.225108e-06
## sd 0.1036970 0.1036339 0.06204795 11.33341 3.180183e-05
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -451.73375 921.46749 956.90532 243.58870 371.000000
## sd 1.732262 28.85099 55.67972 52.14376 39.19745 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.2744141 0.9606934 0.03930664 0.2744141
## sd 0 0.4462730 0.1943472 0.19434724 0.4462730
## significant_negative
## mean 0
## sd 0
##
## $pers_a
## estimate std.error statistic p.value conf.low conf.high
## mean 0.05484896 0.043310646 1.2488281 0.3161042 -0.03031614 0.14001406
## sd 0.03452171 0.003320092 0.7388706 0.2662508 0.03219300 0.03784622
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3528556 0.3409354 0.8092101 30.95031 0.0001543231
## sd 0.1093621 0.1094004 0.0651510 11.88643 0.0025729140
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -452.3210 922.64193 958.07975 244.62057 371.000000
## sd 1.732262 30.1185 58.24095 54.74094 41.33889 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.1984863 0.9831543 0.0168457 0.1984863
## sd 0 0.3989090 0.1287089 0.1287089 0.3989090
## significant_negative
## mean 0
## sd 0
##
## $pers_n
## estimate std.error statistic p.value conf.low conf.high
## mean -0.01945947 0.043371756 -0.4091114 0.4680019 -0.10474473 0.06582579
## sd 0.04490624 0.003125042 0.9702821 0.2636699 0.04887248 0.04147441
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3511661 0.3392152 0.81027906 30.70800 7.639954e-05
## sd 0.1092723 0.1093028 0.06506678 11.78065 1.248112e-03
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -452.8261 923.65224 959.09007 245.25920 371.000000
## sd 1.732262 30.0687 58.13108 54.60971 41.30493 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.08203125 0.4973145 0.5026855 0
## sd 0 0.27444583 0.5000538 0.5000538 0
## significant_negative
## mean 0.08203125
## sd 0.27444583
plot_all_curves(lm_slope_prev_var, 'lm_slope_prev_var')
## $spec_results_o
##
## $spec_results_c
##
## $spec_results_e
##
## $spec_results_a
##
## $spec_results_n
lm_slope_prev_var_ctrl <- lm(slope_prev_var ~ airport_dist +
conservative + male + age + popdens + manufact +
tourism + academics + medinc + healthcare +
temp + slope_prev_var_lag,
data = df_ger_prev_scaled)
summary(lm_slope_prev_var_ctrl)
##
## Call:
## lm(formula = slope_prev_var ~ airport_dist + conservative + male +
## age + popdens + manufact + tourism + academics + medinc +
## healthcare + temp + slope_prev_var_lag, data = df_ger_prev_scaled)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.4605 -0.4684 -0.0894 0.3536 3.8629
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0107049 0.0370087 -0.289 0.772552
## airport_dist -0.0006466 0.0472676 -0.014 0.989094
## conservative 0.0578272 0.0574076 1.007 0.314452
## male -0.1063530 0.0542440 -1.961 0.050679 .
## age -0.3348173 0.0704594 -4.752 2.90e-06 ***
## popdens -0.1144899 0.0591642 -1.935 0.053745 .
## manufact 0.0844884 0.0428457 1.972 0.049371 *
## tourism -0.1367335 0.0426814 -3.204 0.001476 **
## academics -0.0251356 0.0576842 -0.436 0.663279
## medinc 0.2276969 0.0521603 4.365 1.65e-05 ***
## healthcare -0.0223008 0.0487212 -0.458 0.647423
## temp -0.1940683 0.0505286 -3.841 0.000145 ***
## slope_prev_var_lag 0.4329906 0.0602916 7.182 3.88e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7199 on 366 degrees of freedom
## Multiple R-squared: 0.4982, Adjusted R-squared: 0.4817
## F-statistic: 30.28 on 12 and 366 DF, p-value: < 2.2e-16
Predict Cases
lm_cases_0331 <- spec_calculate(df = df_ger_death_scaled,
y = "case_rate_0331",
model = "lm",
controls = covariates %>% append('case_rate_0331_lag'),
all.comb = T)
calc_all_summaries(lm_cases_0331)
## $pers_o
## estimate std.error statistic p.value conf.low conf.high
## mean 0.03158699 0.051022780 0.5994671 0.5627758 -0.06874312 0.13191710
## sd 0.04075752 0.003039436 0.7655479 0.2958935 0.03824516 0.04394422
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.27125522 0.2577158 0.86029794 20.714339 0.0002426819
## sd 0.08294606 0.0823706 0.04660655 7.515557 0.0038974652
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -476.16370 970.32739 1005.76522 275.46553 371.000000
## sd 1.732262 20.54957 39.18622 36.10197 31.35361 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.06103516 0.7888184 0.2111816 0.06103516
## sd 0 0.23942402 0.4081968 0.4081968 0.23942402
## significant_negative
## mean 0
## sd 0
##
## $pers_c
## estimate std.error statistic p.value conf.low conf.high
## mean -0.004091242 0.045901847 -0.08596589 0.5796038 -0.09435164 0.08616916
## sd 0.034700872 0.002269187 0.74910889 0.2613804 0.03536478 0.03460409
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.27028449 0.25673269 0.8607943 20.626725 0.002355802
## sd 0.08560184 0.08506666 0.0479699 7.724115 0.031524022
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -476.35290 970.70581 1006.14364 275.8325 371.000000
## sd 1.732262 21.07044 40.23365 37.14286 32.3575 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0 0.4645996 0.5354004 0
## sd 0 0 0.4988061 0.4988061 0
## significant_negative
## mean 0
## sd 0
##
## $pers_e
## estimate std.error statistic p.value conf.low conf.high
## mean 0.08714459 0.046581747 1.8422270 0.1520497 -0.004452756 0.17874193
## sd 0.04647008 0.001763032 0.9112771 0.1493776 0.043558178 0.04945351
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.27812366 0.26469914 0.85640497 21.418135 2.574047e-06
## sd 0.07660416 0.07592723 0.04326494 7.115188 3.140562e-05
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -474.51363 967.0273 1002.46509 272.86926 371.000000
## sd 1.732262 19.24243 36.5560 33.48926 28.95637 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.3618164 1 0 0.3618164
## sd 0 0.4805847 0 0 0.4805847
## significant_negative
## mean 0
## sd 0
##
## $pers_a
## estimate std.error statistic p.value conf.low conf.high
## mean 0.07275316 0.045946320 1.5732114 0.1733925 -0.01759469 0.16310101
## sd 0.02970329 0.002372876 0.5983738 0.1394271 0.02789450 0.03209363
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.27494180 0.26146891 0.85812315 21.120691 0.0001372159
## sd 0.08246331 0.08190365 0.04643572 7.589597 0.0017659528
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -475.20596 968.41191 1003.8497 274.07200 371.000000
## sd 1.732262 20.51326 39.13057 36.0854 31.17113 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.2634277 1 0 0.2634277
## sd 0 0.4405462 0 0 0.4405462
## significant_negative
## mean 0
## sd 0
##
## $pers_n
## estimate std.error statistic p.value conf.low conf.high
## mean -0.06299523 0.045987096 -1.352350 0.2837433 -0.15342327 0.02743280
## sd 0.03825584 0.002107881 0.770264 0.2261862 0.04073718 0.03608117
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.27430368 0.26081464 0.85855349 21.04466 2.747966e-05
## sd 0.08078661 0.08020112 0.04551672 7.43662 3.526301e-04
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -475.41660 968.83321 1004.27103 274.31321 371.000000
## sd 1.732262 20.13405 38.37129 35.33896 30.53734 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.2299805 0 1 0
## sd 0 0.4208714 0 0 0
## significant_negative
## mean 0.2299805
## sd 0.4208714
plot_all_curves(lm_cases_0331, 'lm_cases_0331')
## $spec_results_o
##
## $spec_results_c
##
## $spec_results_e
##
## $spec_results_a
##
## $spec_results_n
lm_cases_0331_ctrl <- lm(case_rate_0331 ~ airport_dist +
conservative + male + age + popdens + manufact +
tourism + academics + medinc + healthcare +
temp + case_rate_0331_lag,
data = df_ger_death_scaled)
summary(lm_cases_0331_ctrl)
##
## Call:
## lm(formula = case_rate_0331 ~ airport_dist + conservative + male +
## age + popdens + manufact + tourism + academics + medinc +
## healthcare + temp + case_rate_0331_lag, data = df_ger_death_scaled)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.4908 -0.4246 -0.1322 0.1935 5.7140
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.005574 0.041192 -0.135 0.892434
## airport_dist 0.007239 0.052639 0.138 0.890700
## conservative 0.021561 0.064001 0.337 0.736398
## male -0.079280 0.060385 -1.313 0.190032
## age -0.287805 0.078147 -3.683 0.000265 ***
## popdens -0.094136 0.065815 -1.430 0.153484
## manufact 0.045203 0.047719 0.947 0.344117
## tourism -0.094058 0.047518 -1.979 0.048517 *
## academics -0.011510 0.063625 -0.181 0.856548
## medinc 0.184173 0.059260 3.108 0.002032 **
## healthcare -0.036733 0.054269 -0.677 0.498913
## temp -0.180424 0.055771 -3.235 0.001327 **
## case_rate_0331_lag 0.404243 0.071769 5.633 3.54e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8017 on 366 degrees of freedom
## Multiple R-squared: 0.3777, Adjusted R-squared: 0.3573
## F-statistic: 18.51 on 12 and 366 DF, p-value: < 2.2e-16
lm_cases_0930 <- spec_calculate(df = df_ger_death_scaled,
y = "case_rate_0930",
model = "lm",
controls = covariates %>% append('case_rate_0930_lag'),
all.comb = T)
calc_all_summaries(lm_cases_0930)
## $pers_o
## estimate std.error statistic p.value conf.low conf.high
## mean 0.04703511 0.048759556 0.9328083 0.4213025 -0.04884462 0.14291485
## sd 0.04623303 0.003708065 0.8870397 0.3049813 0.04263197 0.05063419
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3311692 0.3188131 0.82304322 28.08371 0.0000322385
## sd 0.1042954 0.1042904 0.06154394 11.21308 0.0004799432
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -458.89825 935.79650 971.23433 252.81805 371.000000
## sd 1.732262 28.15626 54.43081 51.21137 39.42368 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.1430664 0.8413086 0.1586914 0.1430664
## sd 0 0.3501833 0.3654327 0.3654327 0.3501833
## significant_negative
## mean 0
## sd 0
##
## $pers_c
## estimate std.error statistic p.value conf.low conf.high
## mean -0.002610167 0.044018978 -0.02673773 0.5244260 -0.08916812 0.08394778
## sd 0.040295483 0.003040096 0.89565674 0.2664612 0.04373607 0.03749754
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3292839 0.3169004 0.82408293 27.87911 0.0006614882
## sd 0.1074214 0.1074516 0.06314984 11.42535 0.0139028188
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -459.32868 936.65736 972.09518 253.53069 371.000000
## sd 1.732262 28.77586 55.66624 52.42724 40.60527 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.0324707 0.5195312 0.4804688 0
## sd 0 0.1772682 0.4996794 0.4996794 0
## significant_negative
## mean 0.0324707
## sd 0.1772682
##
## $pers_e
## estimate std.error statistic p.value conf.low conf.high
## mean 0.11873043 0.044310503 2.6373935 0.03352240 0.03159922 0.20586163
## sd 0.04761501 0.002634405 0.9060562 0.04271978 0.04339097 0.05201209
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.34275305 0.33059563 0.81616073 29.48154 1.913204e-08
## sd 0.09612161 0.09600844 0.05733096 10.80161 3.024722e-07
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -455.83276 929.6655 965.10334 248.43935 371.000000
## sd 1.732262 26.55416 51.2221 48.02528 36.33397 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.7509766 1 0 0.7509766
## sd 0 0.4325002 0 0 0.4325002
## significant_negative
## mean 0
## sd 0
##
## $pers_a
## estimate std.error statistic p.value conf.low conf.high
## mean 0.001807263 0.044142706 0.02930601 0.6298433 -0.08499398 0.08860851
## sd 0.029309078 0.003273494 0.63816227 0.2379830 0.02826462 0.03165452
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.328552 0.3161559 0.8245256 27.78624 0.0008432176
## sd 0.107648 0.1076773 0.0632661 11.42428 0.0168326235
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -459.52910 937.0582 972.49602 253.80733 371.000000
## sd 1.732262 28.81963 55.7506 52.50393 40.69093 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0 0.5290527 0.4709473 0
## sd 0 0 0.4992162 0.4992162 0
## significant_negative
## mean 0
## sd 0
##
## $pers_n
## estimate std.error statistic p.value conf.low conf.high
## mean -0.001840958 0.044147329 -0.01261312 0.5589962 -0.08865129 0.08496938
## sd 0.037168325 0.003113077 0.79081236 0.2465463 0.04075222 0.03430969
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.3290563 0.3166673 0.82425385 27.83814 0.0001595196
## sd 0.1066349 0.1066540 0.06276337 11.33910 0.0028612546
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -459.41897 936.83794 972.27576 253.61671 371.000000
## sd 1.732262 28.63422 55.37859 52.13274 40.30798 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.01220703 0.5891113 0.4108887 0
## sd 0 0.10982242 0.4920552 0.4920552 0
## significant_negative
## mean 0.01220703
## sd 0.10982242
plot_all_curves(lm_cases_0930, 'lm_cases_0930')
## $spec_results_o
##
## $spec_results_c
##
## $spec_results_e
##
## $spec_results_a
##
## $spec_results_n
lm_cases_0930_ctrl <- lm(case_rate_0930 ~ airport_dist +
conservative + male + age + popdens + manufact +
tourism + academics + medinc + healthcare +
temp + case_rate_0930_lag,
data = df_ger_death_scaled)
summary(lm_cases_0930_ctrl)
##
## Call:
## lm(formula = case_rate_0930 ~ airport_dist + conservative + male +
## age + popdens + manufact + tourism + academics + medinc +
## healthcare + temp + case_rate_0930_lag, data = df_ger_death_scaled)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6305 -0.4216 -0.1215 0.3040 3.5774
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.003935 0.038008 -0.104 0.917604
## airport_dist -0.045205 0.048579 -0.931 0.352702
## conservative 0.102075 0.058827 1.735 0.083551 .
## male -0.034565 0.055737 -0.620 0.535548
## age -0.350952 0.071996 -4.875 1.63e-06 ***
## popdens 0.152294 0.061336 2.483 0.013478 *
## manufact 0.062070 0.044049 1.409 0.159651
## tourism -0.109885 0.043906 -2.503 0.012760 *
## academics -0.217018 0.059681 -3.636 0.000316 ***
## medinc 0.190831 0.056091 3.402 0.000742 ***
## healthcare 0.013822 0.049981 0.277 0.782285
## temp -0.156873 0.051285 -3.059 0.002385 **
## case_rate_0930_lag 0.479824 0.068459 7.009 1.16e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7399 on 366 degrees of freedom
## Multiple R-squared: 0.47, Adjusted R-squared: 0.4526
## F-statistic: 27.05 on 12 and 366 DF, p-value: < 2.2e-16
Predict Deaths
lm_deaths_0331 <- spec_calculate(df = df_ger_death_scaled,
y = "death_rate_0331",
model = "lm",
controls = covariates %>% append('death_rate_0331_lag'),
all.comb = T)
calc_all_summaries(lm_deaths_0331)
## $pers_o
## estimate std.error statistic p.value conf.low conf.high
## mean 0.05926579 0.056344799 1.048407 0.3391156 -0.05152945 0.17006104
## sd 0.02649895 0.002223663 0.465681 0.2207861 0.02562865 0.02803226
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.11230093 0.09564048 0.95084691 6.817236 0.003741489
## sd 0.03249807 0.03024746 0.01580279 1.831876 0.022251556
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -514.578479 1047.15696 1082.59478 335.55025 371.000000
## sd 1.732262 6.852099 11.52894 9.74724 12.28427 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.02636719 0.9821777 0.01782227 0.02636719
## sd 0 0.16024428 0.1323212 0.13232123 0.16024428
## significant_negative
## mean 0
## sd 0
##
## $pers_c
## estimate std.error statistic p.value conf.low conf.high
## mean -0.02176714 0.0507656910 -0.4272184 0.6489583 -0.12159172 0.07805745
## sd 0.02201730 0.0008507798 0.4327864 0.2340294 0.02246288 0.02169193
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.11002111 0.09331891 0.95206101 6.651521 0.00963241
## sd 0.03320389 0.03098238 0.01615476 1.883209 0.05665206
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -515.059950 1048.11990 1083.557726 336.41202 371.000000
## sd 1.732262 6.973763 11.79201 9.998775 12.55107 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0 0.1816406 0.8183594 0
## sd 0 0 0.3855951 0.3855951 0
## significant_negative
## mean 0
## sd 0
##
## $pers_e
## estimate std.error statistic p.value conf.low conf.high
## mean 0.06219816 0.0515241068 1.2102704 0.2768198 -0.03911776 0.16351408
## sd 0.02481655 0.0007406798 0.4899089 0.1810235 0.02544552 0.02425882
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.11323239 0.09658630 0.95035178 6.898352 0.002436617
## sd 0.03213563 0.02994777 0.01566049 1.848994 0.014602120
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -514.381924 1046.76385 1082.201674 335.19815 371.000000
## sd 1.732262 6.787753 11.45944 9.857845 12.14727 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.09594727 1 0 0.09594727
## sd 0 0.29455487 0 0 0.29455487
## significant_negative
## mean 0
## sd 0
##
## $pers_a
## estimate std.error statistic p.value conf.low conf.high
## mean 0.07596961 0.050837168 1.4956105 0.1615796 -0.02399553 0.17593475
## sd 0.01854741 0.001137139 0.3694743 0.1092908 0.01880931 0.01855325
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.11477591 0.09816324 0.9495095 6.996946 0.003652699
## sd 0.03348276 0.03131897 0.0163866 1.927123 0.022385471
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -514.040955 1046.08191 1081.51974 334.61471 371.000000
## sd 1.732262 7.080364 12.03403 10.26139 12.65648 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.1005859 1 0 0.1005859
## sd 0 0.3008164 0 0 0.3008164
## significant_negative
## mean 0
## sd 0
##
## $pers_n
## estimate std.error statistic p.value conf.low conf.high
## mean -0.05389933 0.050903492 -1.0622345 0.3222672 -0.15399488 0.04619623
## sd 0.01914286 0.001035423 0.3829244 0.1659477 0.01845582 0.02001427
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.11220076 0.09553655 0.95090074 6.820286 0.003463388
## sd 0.03252479 0.03033547 0.01585213 1.860572 0.020486709
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -514.599597 1047.19919 1082.637020 335.58811 371.000000
## sd 1.732262 6.860555 11.59847 9.943315 12.29437 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.005859375 0 1 0
## sd 0 0.076331286 0 0 0
## significant_negative
## mean 0.005859375
## sd 0.076331286
plot_all_curves(lm_deaths_0331, 'lm_deaths_0331')
## $spec_results_o
##
## $spec_results_c
##
## $spec_results_e
##
## $spec_results_a
##
## $spec_results_n
lm_deaths_0331_ctrl <- lm(death_rate_0331 ~ airport_dist +
conservative + male + age + popdens + manufact +
tourism + academics + medinc + healthcare +
temp + death_rate_0331_lag,
data = df_ger_death_scaled)
summary(lm_deaths_0331_ctrl)
##
## Call:
## lm(formula = death_rate_0331 ~ airport_dist + conservative +
## male + age + popdens + manufact + tourism + academics + medinc +
## healthcare + temp + death_rate_0331_lag, data = df_ger_death_scaled)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.8867 -0.4742 -0.2131 0.1644 5.1407
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.003276 0.047481 -0.069 0.945030
## airport_dist 0.063879 0.060650 1.053 0.292924
## conservative 0.033146 0.073530 0.451 0.652415
## male -0.170458 0.069611 -2.449 0.014805 *
## age -0.282595 0.087552 -3.228 0.001360 **
## popdens -0.056660 0.075956 -0.746 0.456170
## manufact 0.190338 0.055019 3.460 0.000605 ***
## tourism -0.071972 0.054777 -1.314 0.189699
## academics -0.081969 0.072606 -1.129 0.259654
## medinc 0.061189 0.065733 0.931 0.352534
## healthcare -0.095838 0.062240 -1.540 0.124469
## temp -0.120116 0.064268 -1.869 0.062423 .
## death_rate_0331_lag 0.231548 0.087866 2.635 0.008766 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.924 on 366 degrees of freedom
## Multiple R-squared: 0.1733, Adjusted R-squared: 0.1462
## F-statistic: 6.392 on 12 and 366 DF, p-value: 2.503e-10
lm_deaths_0930 <- spec_calculate(df = df_ger_death_scaled,
y = "death_rate_0930",
model = "lm",
controls = covariates %>% append('death_rate_0930_lag'),
all.comb = T)
calc_all_summaries(lm_deaths_0930)
## $pers_o
## estimate std.error statistic p.value conf.low conf.high
## mean 0.09176990 0.055963743 1.6313826 0.145225 -0.01827605 0.20181584
## sd 0.02910357 0.002231346 0.4896016 0.137522 0.02703758 0.03164673
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.12375372 0.1073184 0.94465296 7.626281 0.002057778
## sd 0.03576864 0.0336473 0.01767784 2.103161 0.013217159
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -512.087716 1042.17543 1077.6133 331.22109 371.000000
## sd 1.732262 7.626127 13.09662 11.0166 13.52055 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.2592773 1 0 0.2592773
## sd 0 0.4382916 0 0 0.4382916
## significant_negative
## mean 0
## sd 0
##
## $pers_c
## estimate std.error statistic p.value conf.low conf.high
## mean -0.02373502 0.0505987501 -0.4672378 0.6381103 -0.12323134 0.07576129
## sd 0.02148019 0.0009822833 0.4226164 0.2372276 0.02203167 0.02109180
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.11782829 0.10128121 0.94783532 7.197791 0.009895088
## sd 0.03662981 0.03456711 0.01808525 2.168291 0.058343172
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -513.359819 1044.71964 1080.15746 333.46091 371.000000
## sd 1.732262 7.745196 13.38859 11.39777 13.84607 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0 0.1352539 0.8647461 0
## sd 0 0 0.3420363 0.3420363 0
## significant_negative
## mean 0
## sd 0
##
## $pers_e
## estimate std.error statistic p.value conf.low conf.high
## mean 0.05256399 0.051225663 1.0284875 0.3487208 -0.04816508 0.15329305
## sd 0.02270246 0.000942625 0.4474427 0.1962162 0.02322695 0.02232001
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.11989229 0.10338102 0.94673451 7.358761 0.003794686
## sd 0.03579107 0.03374005 0.01768661 2.135709 0.021709173
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -512.922157 1043.8443 1079.28214 332.68071 371.000000
## sd 1.732262 7.596441 13.1105 11.22513 13.52902 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.04541016 1 0 0.04541016
## sd 0 0.20822742 0 0 0.20822742
## significant_negative
## mean 0
## sd 0
##
## $pers_a
## estimate std.error statistic p.value conf.low conf.high
## mean 0.01004584 0.05075540 0.1959914 0.7221615 -0.08975852 0.10985019
## sd 0.02082310 0.00116779 0.4083976 0.1840059 0.02043765 0.02144881
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.11736593 0.1008111 0.94808125 7.161593 0.01013347
## sd 0.03683933 0.0347719 0.01819243 2.173253 0.06061994
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -513.457345 1044.91469 1080.3525 333.63568 371.000000
## sd 1.732262 7.789187 13.46958 11.4449 13.92527 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0 0.6926270 0.3073730 0
## sd 0 0 0.4614616 0.4614616 0
## significant_negative
## mean 0
## sd 0
##
## $pers_n
## estimate std.error statistic p.value conf.low conf.high
## mean 0.03733565 0.050718234 0.7335759 0.4914366 -0.06239562 0.13706691
## sd 0.02251156 0.001073727 0.4357044 0.2497344 0.02195929 0.02324317
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.11858409 0.10205446 0.94741952 7.243279 0.01153319
## sd 0.03739883 0.03533075 0.01848994 2.202537 0.06886308
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -513.190383 1044.38077 1079.81859 333.17521 371.000000
## sd 1.732262 7.912561 13.70082 11.59235 14.13676 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0 0.9448242 0.05517578 0
## sd 0 0 0.2283509 0.22835092 0
## significant_negative
## mean 0
## sd 0
plot_all_curves(lm_deaths_0930, 'lm_deaths_0930')
## $spec_results_o
##
## $spec_results_c
##
## $spec_results_e
##
## $spec_results_a
##
## $spec_results_n
lm_deaths_0930_ctrl <- lm(death_rate_0930 ~ airport_dist +
conservative + male + age + popdens + manufact +
tourism + academics + medinc + healthcare +
temp + death_rate_0930_lag,
data = df_ger_death_scaled)
summary(lm_deaths_0930_ctrl)
##
## Call:
## lm(formula = death_rate_0930 ~ airport_dist + conservative +
## male + age + popdens + manufact + tourism + academics + medinc +
## healthcare + temp + death_rate_0930_lag, data = df_ger_death_scaled)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.9553 -0.5092 -0.2248 0.2556 6.5021
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.004634 0.047233 -0.098 0.92190
## airport_dist 0.026542 0.060328 0.440 0.66023
## conservative 0.026796 0.073079 0.367 0.71407
## male -0.087034 0.069386 -1.254 0.21052
## age -0.211738 0.086909 -2.436 0.01531 *
## popdens 0.029944 0.075712 0.396 0.69270
## manufact 0.135993 0.054799 2.482 0.01353 *
## tourism -0.093115 0.054578 -1.706 0.08884 .
## academics -0.162248 0.072481 -2.238 0.02579 *
## medinc 0.093815 0.065408 1.434 0.15234
## healthcare -0.026707 0.061800 -0.432 0.66588
## temp -0.184670 0.064375 -2.869 0.00436 **
## death_rate_0930_lag 0.284429 0.087564 3.248 0.00127 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9191 on 366 degrees of freedom
## Multiple R-squared: 0.182, Adjusted R-squared: 0.1552
## F-statistic: 6.788 on 12 and 366 DF, p-value: 4.504e-11
Predict Socdist Onset
cox_onset_socdist <- spec_calculate(df = df_ger_socdist_scaled,
y = "Surv(cpt_day_socdist_2, event)",
model = "cox_model",
controls = covariates %>%
append(c('cpt_day_socdist_2_lag', 'onset_prev', 'slope_prev')),
all.comb = T)
cox_onset_socdist <- cox_onset_socdist %>% map(filter_socdist)
cox_onset_socdist_hr <- cox_onset_socdist %>% map(coef_to_hr)
calc_all_summaries(cox_onset_socdist_hr)
## $pers_o
## estimate std.error statistic p.value conf.low conf.high fit_n
## mean 0.92412115 0.060146732 -1.3205350 0.2201715 0.82135075 1.03977081 379
## sd 0.02156262 0.002248289 0.4019257 0.1387493 0.01899451 0.02532624 0
## fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean 379 38.32341 0.005706523 36.55378
## sd 0 13.75696 0.024395795 12.35949
## fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean 0.005739844 36.14392 0.005735019 NA
## sd 0.023985584 12.05654 0.023777053 NA
## fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean NA 0.09557573 0.9999496 0.61749673
## sd NA 0.03296967 0.0000000 0.03336312
## fit_std.error.concordance fit_logLik fit_AIC fit_BIC fit_nobs
## mean 0.028868295 -1856.052442 3730.10488 3765.543 379
## sd 0.001210389 6.878482 12.02107 11.250 0
## significant positive negative significant_positive significant_negative
## mean 0.06811523 0 1 0 0.06811523
## sd 0.25197430 0 0 0 0.25197430
##
## $pers_c
## estimate std.error statistic p.value conf.low conf.high fit_n
## mean 0.92379825 0.053425612 -1.4965299 0.1765796 0.8319259 1.02582137 379
## sd 0.02216087 0.001129165 0.4722299 0.1393372 0.0186606 0.02630592 0
## fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean 379 38.89125 0.002080154 37.09716
## sd 0 12.54155 0.007141225 11.33829
## fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean 0.002219141 36.66871 0.002244855 NA
## sd 0.007425369 11.03882 0.007438921 NA
## fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean NA 0.09703073 0.9999496 0.62825725
## sd NA 0.02998817 0.0000000 0.02806221
## fit_std.error.concordance fit_logLik fit_AIC fit_BIC fit_nobs
## mean 0.02888723 -1855.768524 3729.53705 3764.97487 379
## sd 0.00105630 6.270773 10.78844 10.29674 0
## significant positive negative significant_positive significant_negative
## mean 0.1914062 0 1 0 0.1914062
## sd 0.3934561 0 0 0 0.3934561
##
## $pers_e
## estimate std.error statistic p.value conf.low conf.high fit_n
## mean 0.88750828 0.0562011720 -2.1242714 0.04010842 0.79494811 0.99084875 379
## sd 0.01407696 0.0008904741 0.2702866 0.02493297 0.01325974 0.01506055 0
## fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean 379 41.01862 0.002206119 39.25972
## sd 0 13.53017 0.009699107 12.25642
## fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean 0.00234413 38.82729 0.002322961 NA
## sd 0.01021985 11.92044 0.010047201 NA
## fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean NA 0.1020040 0.9999496 0.62851368
## sd NA 0.0321979 0.0000000 0.02755481
## fit_std.error.concordance fit_logLik fit_AIC fit_BIC fit_nobs
## mean 0.028252338 -1854.704839 3727.40968 3762.84750 379
## sd 0.001407939 6.765084 11.72497 10.85692 0
## significant positive negative significant_positive significant_negative
## mean 0.6994629 0 1 0 0.6994629
## sd 0.4585476 0 0 0 0.4585476
##
## $pers_a
## estimate std.error statistic p.value conf.low conf.high fit_n
## mean 0.94998162 0.055857806 -0.9252410 0.3798245 0.85145243 1.0599184 379
## sd 0.01686576 0.001206494 0.3290921 0.1672532 0.01424408 0.0201586 0
## fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean 379 37.40628 0.00487520 35.73613
## sd 0 13.19216 0.01728311 11.95547
## fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean 0.005074892 35.31908 0.005112128 NA
## sd 0.017780079 11.64557 0.017763039 NA
## fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean NA 0.09343289 0.9999496 0.62513121
## sd NA 0.03167007 0.0000000 0.03199302
## fit_std.error.concordance fit_logLik fit_AIC fit_BIC fit_nobs
## mean 0.028872705 -1856.51101 3731.02201 3766.45984 379
## sd 0.001126387 6.59608 11.46145 10.84048 0
## significant positive negative significant_positive significant_negative
## mean 0 0 1 0 0
## sd 0 0 0 0 0
##
## $pers_n
## estimate std.error statistic p.value conf.low conf.high fit_n
## mean 1.08916744 0.057407948 1.4897329 0.14528486 0.97327553 1.21886580 379
## sd 0.01220407 0.001200459 0.2138959 0.06005665 0.01235991 0.01223976 0
## fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean 379 38.68176 0.00353503 37.02830
## sd 0 13.43207 0.01343931 12.20186
## fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean 0.003639324 36.63856 0.003605216 NA
## sd 0.013488312 11.86728 0.013282388 NA
## fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean NA 0.09645849 0.9999496 0.62898123
## sd NA 0.03214506 0.0000000 0.03045481
## fit_std.error.concordance fit_logLik fit_AIC fit_BIC fit_nobs
## mean 0.028670041 -1855.873268 3729.74654 3765.18436 379
## sd 0.001247607 6.716033 11.72373 11.09026 0
## significant positive negative significant_positive significant_negative
## mean 0 1 0 0 0
## sd 0 0 0 0 0
plot_all_curves(cox_onset_socdist_hr, 'cox_onset_socdist_hr', hr = T)
## $spec_results_o
##
## $spec_results_c
##
## $spec_results_e
##
## $spec_results_a
##
## $spec_results_n
cox_onset_socdist_ctrl <- coxph(Surv(cpt_day_socdist_2, event) ~ airport_dist +
conservative + male + age + popdens + manufact +
tourism + academics + medinc + healthcare +
temp + onset_prev + slope_prev +
cpt_day_socdist_2_lag,
data = df_ger_socdist_scaled)
summary(cox_onset_socdist_ctrl)
## Call:
## coxph(formula = Surv(cpt_day_socdist_2, event) ~ airport_dist +
## conservative + male + age + popdens + manufact + tourism +
## academics + medinc + healthcare + temp + onset_prev + slope_prev +
## cpt_day_socdist_2_lag, data = df_ger_socdist_scaled)
##
## n= 379, number of events= 379
##
## coef exp(coef) se(coef) z Pr(>|z|)
## airport_dist 0.125877 1.134142 0.067907 1.854 0.0638 .
## conservative -0.331611 0.717767 0.082851 -4.002 6.27e-05 ***
## male -0.141864 0.867740 0.081693 -1.737 0.0825 .
## age -0.017427 0.982724 0.103089 -0.169 0.8658
## popdens 0.108749 1.114883 0.091485 1.189 0.2346
## manufact -0.072176 0.930367 0.052325 -1.379 0.1678
## tourism 0.009962 1.010012 0.061504 0.162 0.8713
## academics -0.058010 0.943640 0.079147 -0.733 0.4636
## medinc -0.104273 0.900979 0.071349 -1.461 0.1439
## healthcare -0.173877 0.840400 0.069290 -2.509 0.0121 *
## temp 0.077350 1.080420 0.072919 1.061 0.2888
## onset_prev -0.013565 0.986527 0.075247 -0.180 0.8569
## slope_prev 0.146099 1.157311 0.086200 1.695 0.0901 .
## cpt_day_socdist_2_lag -0.238041 0.788170 0.114252 -2.083 0.0372 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## airport_dist 1.1341 0.8817 0.9928 1.2956
## conservative 0.7178 1.3932 0.6102 0.8443
## male 0.8677 1.1524 0.7394 1.0184
## age 0.9827 1.0176 0.8029 1.2028
## popdens 1.1149 0.8970 0.9319 1.3338
## manufact 0.9304 1.0748 0.8397 1.0308
## tourism 1.0100 0.9901 0.8953 1.1394
## academics 0.9436 1.0597 0.8080 1.1020
## medinc 0.9010 1.1099 0.7834 1.0362
## healthcare 0.8404 1.1899 0.7337 0.9626
## temp 1.0804 0.9256 0.9365 1.2464
## onset_prev 0.9865 1.0137 0.8513 1.1433
## slope_prev 1.1573 0.8641 0.9774 1.3703
## cpt_day_socdist_2_lag 0.7882 1.2688 0.6300 0.9860
##
## Concordance= 0.678 (se = 0.029 )
## Likelihood ratio test= 61.37 on 14 df, p=7e-08
## Wald test = 57.4 on 14 df, p=3e-07
## Score (logrank) test = 58.42 on 14 df, p=2e-07
Predict Socdist Mean
lm_mean_socdist <- spec_calculate(df = df_ger_socdist_scaled,
y = "mean_diff_socdist_2",
model = "lm",
controls = covariates %>%
append(c('mean_diff_socdist_2_lag', 'onset_prev', 'slope_prev')),
all.comb = T)
lm_mean_socdist <- lm_mean_socdist %>% map(filter_socdist)
calc_all_summaries(lm_mean_socdist)
## $pers_o
## estimate std.error statistic p.value conf.low conf.high
## mean 0.15330885 0.044355910 3.425476 0.02297009 0.06608679 0.24053090
## sd 0.07755911 0.001414087 1.656772 0.03632131 0.07585609 0.07932302
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.45206387 0.43888864 0.7476768 35.209876 7.055131e-20
## sd 0.07030338 0.07031662 0.0457296 7.949236 2.044565e-18
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 9.000000 -421.8148 865.62962 908.9425 207.11986 369.000000
## sd 1.732262 23.1611 44.31669 40.9420 26.57468 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.8122559 1 0 0.8122559
## sd 0 0.3905554 0 0 0.3905554
## significant_negative
## mean 0
## sd 0
##
## $pers_c
## estimate std.error statistic p.value conf.low conf.high
## mean -0.03851258 0.04069143 -0.9421542 0.4218812 -0.11852873 0.04150356
## sd 0.02790481 0.00308532 0.6560585 0.2885431 0.02975273 0.02730813
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.43015918 0.41651193 0.76145621 32.660794 8.760185e-07
## sd 0.09616314 0.09669053 0.06060868 9.707367 2.157079e-05
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 9.000000 -428.28432 878.56865 921.88154 215.39983 369.000000
## sd 1.732262 29.50691 57.06175 53.66655 36.34967 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.0769043 0.0546875 0.9453125 0
## sd 0 0.2664721 0.2273970 0.2273970 0
## significant_negative
## mean 0.0769043
## sd 0.2664721
##
## $pers_e
## estimate std.error statistic p.value conf.low conf.high
## mean -0.01503002 0.041886341 -0.3896659 0.5811354 -0.09739585 0.06733580
## sd 0.02278655 0.003250974 0.5237195 0.1998854 0.01851793 0.02787933
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.42867767 0.41499876 0.76240500 32.47312 3.773528e-06
## sd 0.09720623 0.09774027 0.06116183 9.72353 8.896818e-05
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 9.000000 -428.73869 879.4774 922.79028 215.95984 369.000000
## sd 1.732262 29.73164 57.4969 54.06717 36.74395 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0 0.1518555 0.8481445 0
## sd 0 0 0.3589246 0.3589246 0
## significant_negative
## mean 0
## sd 0
##
## $pers_a
## estimate std.error statistic p.value conf.low conf.high
## mean -0.02332656 0.040841977 -0.5897552 0.4122219 -0.10363874 0.05698563
## sd 0.04379788 0.002890095 1.0885441 0.2944594 0.04279858 0.04549033
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.43026173 0.41662073 0.76130558 32.712822 6.230705e-06
## sd 0.09788511 0.09844834 0.06159559 9.858796 1.500141e-04
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 9.000000 -428.17439 878.3488 921.66168 215.36107 369.000000
## sd 1.732262 29.92968 57.9119 54.51971 37.00057 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0.1582031 0.3818359 0.6181641 0
## sd 0 0.3649759 0.4858960 0.4858960 0
## significant_negative
## mean 0.1582031
## sd 0.3649759
##
## $pers_n
## estimate std.error statistic p.value conf.low conf.high
## mean -0.00604797 0.041254131 -0.1365582 0.7023933 -0.08717062 0.07507468
## sd 0.01968351 0.003240101 0.4720721 0.2018550 0.02250371 0.01869854
## fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean 0.42841045 0.41472473 0.76258917 32.435635 6.023785e-06
## sd 0.09718258 0.09771424 0.06110589 9.710442 1.448236e-04
## fit_df fit_logLik fit_AIC fit_BIC fit_deviance fit_df.residual
## mean 9.000000 -428.83354 879.66708 922.97998 216.06085 369.000000
## sd 1.732262 29.68545 57.40582 53.97958 36.73502 1.732262
## fit_nobs significant positive negative significant_positive
## mean 379 0 0.4108887 0.5891113 0
## sd 0 0 0.4920552 0.4920552 0
## significant_negative
## mean 0
## sd 0
plot_all_curves(lm_mean_socdist, 'lm_mean_socdist')
## $spec_results_o
##
## $spec_results_c
##
## $spec_results_e
##
## $spec_results_a
##
## $spec_results_n
lm_mean_socdist_ctrl <- lm(mean_diff_socdist_2 ~ airport_dist +
conservative + male + age + popdens + manufact +
tourism + academics + medinc + healthcare +
temp + onset_prev + slope_prev +
mean_diff_socdist_2_lag,
data = df_ger_socdist_scaled)
summary(lm_mean_socdist_ctrl)
##
## Call:
## lm(formula = mean_diff_socdist_2 ~ airport_dist + conservative +
## male + age + popdens + manufact + tourism + academics + medinc +
## healthcare + temp + onset_prev + slope_prev + mean_diff_socdist_2_lag,
## data = df_ger_socdist_scaled)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.85241 -0.44592 0.01798 0.37505 2.81498
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.001481 0.035208 -0.042 0.9665
## airport_dist -0.010038 0.045187 -0.222 0.8243
## conservative -0.273373 0.055239 -4.949 1.14e-06 ***
## male -0.117690 0.051800 -2.272 0.0237 *
## age 0.003341 0.067934 0.049 0.9608
## popdens 0.328194 0.056804 5.778 1.63e-08 ***
## manufact 0.085452 0.041025 2.083 0.0380 *
## tourism -0.013999 0.041163 -0.340 0.7340
## academics 0.364853 0.054490 6.696 8.15e-11 ***
## medinc -0.073709 0.050032 -1.473 0.1416
## healthcare -0.008339 0.046183 -0.181 0.8568
## temp -0.081000 0.048476 -1.671 0.0956 .
## onset_prev 0.008089 0.051930 0.156 0.8763
## slope_prev 0.031435 0.055945 0.562 0.5745
## mean_diff_socdist_2_lag 0.135695 0.064960 2.089 0.0374 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6853 on 364 degrees of freedom
## Multiple R-squared: 0.5478, Adjusted R-squared: 0.5304
## F-statistic: 31.49 on 14 and 364 DF, p-value: < 2.2e-16
Save all results
results_ger <- list(cox_onset_prev_hr, lm_slope_prev, lm_slope_prev_var,
lm_cases_0331, lm_cases_0930, lm_deaths_0331, lm_deaths_0930,
cox_onset_socdist_hr, lm_mean_socdist)
names(results_ger) <- list('cox_onset_prev_hr', 'lm_slope_prev', 'lm_slope_prev_var',
'lm_cases_0331', 'lm_cases_0930', 'lm_deaths_0331', 'lm_deaths_0930',
'cox_onset_socdist_hr', 'lm_mean_socdist')
save(results_ger, file="GER_spec_results.RData")
Descriptives
Distributions
df_ger_desc <- df_ger_socdist %>%
dplyr::select(pers_o, pers_c, pers_e, pers_a, pers_n,
age, male, conservative,
academics, medinc, manufact,
airport_dist, tourism, healthcare, popdens) %>%
distinct() %>% drop_na()
ger_means <- df_ger_desc %>% summarise_all(mean)
ger_sd <- df_ger_desc %>% summarise_all(sd)
desc_ger <- rbind(ger_means, ger_sd) %>% t() %>% round(3) %>% as.data.frame()
desc_ger %>% rownames_to_column() %>% write_csv('ger_descriptives.csv')
desc_ger
Explore correlations
a <- df_ger_socdist_unscaled %>% select(kreis, pers_o:pers_n, cpt_day_socdist, mean_diff_socdist)
b <- df_ger_prev_unscaled %>% select(kreis, onset_prev, slope_prev)
c <- df_ger_death_unscaled %>% select(kreis, case_rate_0331:death_rate_0930)
ger_joined <- plyr::join_all(list(a, b, c), by='kreis')
ger_joined %>% select(-kreis) %>% cor(use = 'pairwise.complete')
## pers_o pers_c pers_e pers_a pers_n
## pers_o 1.00000000 -0.10201270 0.27983871 0.15759347 -0.0640347187
## pers_c -0.10201270 1.00000000 0.21472242 0.32690021 -0.3551086811
## pers_e 0.27983871 0.21472242 1.00000000 0.18288787 -0.3649780097
## pers_a 0.15759347 0.32690021 0.18288787 1.00000000 -0.4280315714
## pers_n -0.06403472 -0.35510868 -0.36497801 -0.42803157 1.0000000000
## cpt_day_socdist -0.10416238 0.10437674 -0.04952867 -0.01672102 0.0138442293
## mean_diff_socdist 0.31429946 -0.12682312 0.09840937 -0.08341683 0.0006497013
## onset_prev -0.16574985 0.05702967 -0.27615283 -0.09470167 0.2352513229
## slope_prev 0.15470785 -0.04204164 0.24306903 0.13157797 -0.1983834025
## case_rate_0331 0.10695243 -0.01752091 0.20572870 0.13949525 -0.1722746724
## death_rate_0331 0.06566400 -0.03556727 0.10805800 0.09283184 -0.0999658180
## case_rate_0930 0.15634007 -0.05283919 0.25572247 0.04675307 -0.1116940920
## death_rate_0930 0.07587651 -0.03740828 0.08958204 0.02599047 -0.0095583322
## cpt_day_socdist mean_diff_socdist onset_prev slope_prev
## pers_o -0.1041623824 0.3142994608 -0.16574985 0.15470785
## pers_c 0.1043767434 -0.1268231184 0.05702967 -0.04204164
## pers_e -0.0495286726 0.0984093720 -0.27615283 0.24306903
## pers_a -0.0167210202 -0.0834168254 -0.09470167 0.13157797
## pers_n 0.0138442293 0.0006497013 0.23525132 -0.19838340
## cpt_day_socdist 1.0000000000 0.0358262655 0.14392128 -0.11996515
## mean_diff_socdist 0.0358262655 1.0000000000 -0.20519894 0.22948985
## onset_prev 0.1439212762 -0.2051989385 1.00000000 -0.67617865
## slope_prev -0.1199651501 0.2294898547 -0.67617865 1.00000000
## case_rate_0331 -0.0689721829 0.1733238658 -0.66024501 0.92767374
## death_rate_0331 0.0003684063 0.1435286041 -0.40795429 0.69587384
## case_rate_0930 -0.1260015001 0.2185131719 -0.58477976 0.82060844
## death_rate_0930 -0.0078242968 0.1456354879 -0.32697005 0.64249969
## case_rate_0331 death_rate_0331 case_rate_0930 death_rate_0930
## pers_o 0.10695243 0.0656640020 0.15634007 0.075876508
## pers_c -0.01752091 -0.0355672682 -0.05283919 -0.037408282
## pers_e 0.20572870 0.1080579960 0.25572247 0.089582044
## pers_a 0.13949525 0.0928318364 0.04675307 0.025990472
## pers_n -0.17227467 -0.0999658180 -0.11169409 -0.009558332
## cpt_day_socdist -0.06897218 0.0003684063 -0.12600150 -0.007824297
## mean_diff_socdist 0.17332387 0.1435286041 0.21851317 0.145635488
## onset_prev -0.66024501 -0.4079542869 -0.58477976 -0.326970047
## slope_prev 0.92767374 0.6958738421 0.82060844 0.642499692
## case_rate_0331 1.00000000 0.7967609328 0.82553862 0.687015730
## death_rate_0331 0.79676093 1.0000000000 0.68600498 0.850162063
## case_rate_0930 0.82553862 0.6860049833 1.00000000 0.745023167
## death_rate_0930 0.68701573 0.8501620632 0.74502317 1.000000000
ger_means <- ger_joined %>% select(-kreis) %>% summarise_all(mean)
ger_sd <- ger_joined %>% select(-kreis) %>% summarise_all(sd)
desc_ger <- rbind(ger_means, ger_sd) %>% t() %>% round(3) %>% as.data.frame() %>%
rename(mean=V1, sd=V2)
desc_ger %>% rownames_to_column() %>% write_csv('ger_descriptives.csv')
desc_ger
Visualizations
df_ger_prev <- read_csv('GER_prevalence.csv')
df_ger_prev <- df_ger_prev %>%
mutate(date = as.Date(date, "%d%b%Y"),
kreis = as.character(kreis)) %>%
dplyr::select(kreis, date, rate_day) %>%
filter(date < '2020-04-01')
# create sequence of dates
date_sequence <- seq.Date(min(df_ger_prev$date),
max(df_ger_prev$date), 1)
# create data frame with time sequence
df_dates = data.frame(date_sequence, 1:length(date_sequence))
names(df_dates) <- c('date', 'time')
tto <- df_ger_prev_unscaled %>%
filter(event==1) %>%
merge(df_dates, by.x = 'onset_prev', by.y = 'time') %>%
ggplot(aes(x=date+1)) +
geom_histogram(bins=15) +
theme_light() +
ggtitle('Distribution - COVID-19 Onset') +
theme(plot.title = element_text(size=11)) +
ylab('Absolute frequency') +
xlab('Date')
soa <- df_ger_slope_prev %>%
ggplot(aes(x=slope_prev)) +
geom_histogram(bins=15) +
theme_light() +
ggtitle('Distribution - COVID-19 Growth Rates') +
theme(plot.title = element_text(size=11)) +
ylab('Absolute frequency') +
xlab('Standardized growth rates')
figure <- ggarrange(tto, soa,
labels = c("", ""),
ncol = 2, nrow = 1)
figure
tto <- df_ger_death %>%
ggplot(aes(x=case_rate_0930)) +
geom_histogram(bins=15) +
theme_light() +
ggtitle('Distribution - Cumulative Case Rates') +
theme(plot.title = element_text(size=11)) +
ylab('Absolute frequency') +
xlab('Cases Per 1000 Inhabitants')
soa <- df_ger_death %>%
ggplot(aes(x=death_rate_0930)) +
geom_histogram(bins=15) +
theme_light() +
ggtitle('Distribution - Cumulative Death Rates') +
theme(plot.title = element_text(size=11)) +
ylab('Absolute frequency') +
xlab('Deaths Per 1000 Inhabitants')
figure <- ggarrange(tto, soa,
labels = c("", ""),
ncol = 2, nrow = 1)
figure
# create sequence of dates
date_sequence <- seq.Date(min(df_ger_socdist$date),
max(df_ger_socdist$date), 1)
# create data frame with time sequence
df_dates = data.frame(date_sequence, 1:length(date_sequence))
names(df_dates) <- c('date', 'time')
tto <- df_ger_cpt_socdist %>%
merge(df_dates, by.x = 'cpt_day_socdist_2', by.y = 'time') %>%
#filter(cpt_day_socdist_2 < 30) %>%
ggplot(aes(x=date+1)) +
geom_histogram(bins=15) +
theme_light() +
ggtitle('Distribution - Time to Adoption') +
theme(plot.title = element_text(size=11)) +
ylab('Absolute frequency') +
xlab('Date') +
xlim(c(as.Date('2020-03-01'), as.Date('2020-03-31')))
soa <- df_ger_cpt_socdist %>%
ggplot(aes(x=mean_diff_socdist_2)) +
geom_histogram(bins=15) +
theme_light() +
ggtitle('Distribution - Strength of Adjustment') +
theme(plot.title = element_text(size=11)) +
ylab('Absolute frequency') +
xlab('Standardized mean difference')
figure <- ggarrange(tto, soa,
labels = c("", ""),
ncol = 2, nrow = 1)
figure
>>>>>>> f9bb205265c8cd2a54b17b2b2634fb89e3d3c38f
Social Distancing - Change point analysis